Research Script
Eze-Emiri CN
2025-01-11
Install packages and load libraries
Install required packages and open their respective libraries: sf, rgdal, and readxl
library(sf)## Warning: package 'sf' was built under R version 4.2.3## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUElibrary(rgdal)## Warning: package 'rgdal' was built under R version 4.2.3## Loading required package: sp## Warning: package 'sp' was built under R version 4.2.3## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
## (status 2 uses the sf package in place of rgdal)## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/xyaze/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Users/xyaze/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:2.0-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.library(readxl)## Warning: package 'readxl' was built under R version 4.2.3library(psych)## Warning: package 'psych' was built under R version 4.2.3library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.2.3##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alphalibrary(Hmisc)## Warning: package 'Hmisc' was built under R version 4.2.3##
## Attaching package: 'Hmisc'## The following object is masked from 'package:psych':
##
## describe## The following objects are masked from 'package:base':
##
## format.pval, unitslibrary(corrplot)## Warning: package 'corrplot' was built under R version 4.2.3## corrplot 0.92 loadedlibrary(spdep)## Warning: package 'spdep' was built under R version 4.2.3## Loading required package: spData## Warning: package 'spData' was built under R version 4.2.3library(GWmodel)## Warning: package 'GWmodel' was built under R version 4.2.3## Loading required package: robustbase## Warning: package 'robustbase' was built under R version 4.2.3## Loading required package: Rcpp## Warning: package 'Rcpp' was built under R version 4.2.3## Welcome to GWmodel version 2.3-1.library(tmap)## Warning: package 'tmap' was built under R version 4.2.3## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')library(viridis)## Warning: package 'viridis' was built under R version 4.2.3## Loading required package: viridisLite## Warning: package 'viridisLite' was built under R version 4.2.3library(RColorBrewer)
library(dplyr)## Warning: package 'dplyr' was built under R version 4.2.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:Hmisc':
##
## src, summarize## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(car)## Loading required package: carData## Warning: package 'carData' was built under R version 4.2.3##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:psych':
##
## logitRead Nigeria shape file
Read shape file using readOGR() and rename to ngamap NB: Add geopolitical zone map
ngamap1 <- readOGR(dsn = "C:/Users/xyaze/Documents/NGA_adm", layer = "NGA_adm1", verbose = TRUE)## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\xyaze\Documents\NGA_adm", layer: "NGA_adm1"
## with 38 features
## It has 9 fields
## Integer64 fields read as strings: ID_0 ID_1class(ngamap1)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"head(ngamap1@data) #view columns## ID_0 ISO NAME_0 ID_1 NAME_1 TYPE_1 ENGTYPE_1 NL_NAME_1 VARNAME_1
## 0 163 NGA Nigeria 1 Abia State State <NA> <NA>
## 1 163 NGA Nigeria 2 Adamawa State State <NA> <NA>
## 2 163 NGA Nigeria 3 Akwa Ibom State State <NA> <NA>
## 3 163 NGA Nigeria 4 Anambra State State <NA> <NA>
## 4 163 NGA Nigeria 5 Bauchi State State <NA> <NA>
## 5 163 NGA Nigeria 6 Bayelsa State State <NA> <NA>plot(ngamap1) #visualise the Map of Nigeria
geozone <- read_excel("C:/Users/xyaze/Desktop/Geo_Zone.xlsx", sheet = "Sheet1")
ngageo <- sp::merge(ngamap1, geozone, by = 'NAME_1') #merge two data sets
head(ngageo)## class : SpatialPolygonsDataFrame
## features : 6
## extent : 5.376527, 13.72793, 4.270418, 12.5025 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 10
## names : NAME_1, ID_0, ISO, NAME_0, ID_1, TYPE_1, ENGTYPE_1, NL_NAME_1, VARNAME_1, GEO_ZONE
## min values : Abia, 163, NGA, Nigeria, 1, State, State, NA, NA, North-East
## max values : Bayelsa, 163, NGA, Nigeria, 6, State, State, NA, NA, South-SouthView(ngageo) #view data set in a window
ngageomap <- st_as_sf(ngageo)
class(ngageomap)## [1] "sf" "data.frame"View(ngageomap)
print(ggplot(ngageomap) + geom_sf(aes(fill = GEO_ZONE, color='')) + scale_fill_viridis_d(option = "viridis") + theme_bw()) #plot map
tm_shape(ngageomap) + tm_fill (col = "GEO_ZONE", palette = "viridis", title = "Geopolitical Zones")
##Extract south south data and visualise Plot the map of South-South Nigeria showing the 123 Local Government Areas using the subset() function
ssmap1 <- subset(ngamap1, ID_1 == 3 | ID_1 == 6 | ID_1 == 9 | ID_1 == 10 | ID_1 == 12 | ID_1 == 33) #filter south-south states
ssmap1 <- st_as_sf(ssmap1)
class(ssmap1) #view object type## [1] "sf" "data.frame"head(ssmap1) #snapshot of data -first 6 rows and columns## Simple feature collection with 6 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 4.978159 ymin: 4.270418 xmax: 9.467367 ymax: 7.573423
## Geodetic CRS: WGS 84
## ID_0 ISO NAME_0 ID_1 NAME_1 TYPE_1 ENGTYPE_1 NL_NAME_1 VARNAME_1
## 2 163 NGA Nigeria 3 Akwa Ibom State State <NA> <NA>
## 5 163 NGA Nigeria 6 Bayelsa State State <NA> <NA>
## 8 163 NGA Nigeria 9 Cross River State State <NA> <NA>
## 9 163 NGA Nigeria 10 Delta State State <NA> <NA>
## 11 163 NGA Nigeria 12 Edo State State <NA> <NA>
## 32 163 NGA Nigeria 33 Rivers State State <NA> <NA>
## geometry
## 2 MULTIPOLYGON (((7.610695 4....
## 5 MULTIPOLYGON (((6.440416 4....
## 8 MULTIPOLYGON (((8.981203 6....
## 9 MULTIPOLYGON (((5.391528 5....
## 11 MULTIPOLYGON (((6.029233 7....
## 32 MULTIPOLYGON (((6.970418 4....ggplot(ssmap1) + geom_sf(aes(fill = NAME_1, color='')) + scale_fill_viridis_d(option = "viridis") + theme_bw() #plot map
ngamap <- readOGR(dsn = "C:/Users/xyaze/Documents/NGA_adm", layer = "NGA_adm2", verbose = TRUE)## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\xyaze\Documents\NGA_adm", layer: "NGA_adm2"
## with 775 features
## It has 11 fields
## Integer64 fields read as strings: ID_0 ID_1 ID_2ssmap <- subset(ngamap, ID_1 == 3 | ID_1 == 6 | ID_1 == 9 | ID_1 == 10 | ID_1 == 12 | ID_1 == 33) #filter south-south states
ssmap <- st_as_sf(ssmap)
ggplot(ssmap) + geom_sf(aes(fill = NAME_1, color='')) + scale_fill_viridis_d(option = "viridis") + theme_bw()
tm_shape(ssmap) + tm_fill (col = "NAME_1", palette = "viridis", title = "South-South States")
rschdata <- read_excel("C:/Users/xyaze/Desktop/rschdata.xlsx", sheet = "Sheet1") #import data
head(rschdata) #snapshot of data## # A tibble: 6 × 26
## OBJECTID NAME_1 NAME_2 LGA `AREA (km2)` POP POP_DEN PUBLIC_HF PHF_COV
## <dbl> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 696 Rivers Obio/… Obio… 263. 708008 2691. 27 3.81
## 2 695 Rivers Khana Khana 524. 361163 689. 32 8.86
## 3 69 Akwa Ibom Uyo Uyo 154. 408011 2648. 23 5.64
## 4 41 Akwa Ibom Eket Eket 210. 314900 1502. 19 6.03
## 5 176 Cross Riv… Calab… Cala… 48.8 190551 3905. 73 38.3
## 6 118 Bayelsa Yeneg… Yena… 707. 776979 1099. 45 5.79
## # ℹ 17 more variables: HIV_PREV <dbl>, EST_PLHIV <dbl>, TX_CURR <dbl>,
## # NOR_TXCURR <dbl>, ART_HF <dbl>, AHF_COV <dbl>, `MPI INDEX` <dbl>,
## # DEPRIVATION <dbl>, POP_POV <dbl>, `HEALTH_DEP%` <dbl>, HEALTH_DEP <dbl>,
## # `EDUC_DEP%` <dbl>, EDUC_DEP <dbl>, `WORK_DEP%` <dbl>, WORK_DEP <dbl>,
## # `SEC_SHOCK%` <dbl>, SEC_SHOCK <dbl>Merge spatial and non-spatial data
Merge data sets using similar columns
ssdata <- merge(ssmap, rschdata) #merge two data sets
head(ssdata)## Simple feature collection with 6 features and 35 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 7.465471 ymin: 4.469585 xmax: 8.142541 ymax: 5.164662
## Geodetic CRS: WGS 84
## NAME_1 NAME_2 ID_0 ISO NAME_0 ID_1 ID_2 TYPE_2
## 1 Akwa Ibom Abak 163 NGA Nigeria 3 39 Local Authority
## 2 Akwa Ibom Eastern Obolo 163 NGA Nigeria 3 40 Local Authority
## 3 Akwa Ibom Eket 163 NGA Nigeria 3 41 Local Authority
## 4 Akwa Ibom Esit Eket 163 NGA Nigeria 3 42 Local Authority
## 5 Akwa Ibom Essien-U 163 NGA Nigeria 3 43 Local Authority
## 6 Akwa Ibom EtimEkpo 163 NGA Nigeria 3 44 Local Authority
## ENGTYPE_2 NL_NAME_2 VARNAME_2 OBJECTID LGA AREA (km2)
## 1 Local Authority <NA> <NA> 39 Abak 167.0
## 2 Local Authority <NA> <NA> 40 Eastern Obolo 156.5
## 3 Local Authority <NA> <NA> 41 Eket 209.7
## 4 Local Authority <NA> <NA> 42 Esit - Eket 148.2
## 5 Local Authority <NA> Essien Udim 43 Essien Udim 298.1
## 6 Local Authority <NA> Etim Ekpo 44 Etim Ekpo 200.4
## POP POP_DEN PUBLIC_HF PHF_COV HIV_PREV EST_PLHIV TX_CURR NOR_TXCURR
## 1 154945 927.8144 16 10.326245 5.5 8521.975 5363 3461.228
## 2 125525 802.0767 12 9.559849 5.5 6903.875 5499 4380.801
## 3 314900 1501.6691 19 6.033661 5.5 17319.500 17231 5471.896
## 4 93225 629.0486 15 16.090105 5.5 5127.375 6646 7128.989
## 5 249990 838.6112 21 8.400336 5.5 13749.450 9543 3817.353
## 6 138392 690.5788 14 10.116192 5.5 7611.560 4469 3229.233
## ART_HF AHF_COV MPI INDEX DEPRIVATION POP_POV HEALTH_DEP% HEALTH_DEP
## 1 4 7.458512 0.288 71.3 110475.8 42.0 65076.90
## 2 0 0.000000 0.308 73.2 91884.3 39.4 49456.85
## 3 17 9.865939 0.308 73.2 230506.8 39.4 124070.60
## 4 0 0.000000 0.308 73.2 68240.7 39.4 36730.65
## 5 12 12.574662 0.288 71.3 178242.9 42.0 104995.80
## 6 0 0.000000 0.288 71.3 98673.5 42.0 58124.64
## EDUC_DEP% EDUC_DEP WORK_DEP% WORK_DEP SEC_SHOCK% SEC_SHOCK
## 1 7.7 11930.77 16.4 25410.98 8.0 12395.600
## 2 6.4 8033.60 16.6 20837.15 15.2 19079.800
## 3 6.4 20153.60 20.0 62980.00 5.6 17634.400
## 4 6.4 5966.40 17.4 16221.15 5.4 5034.150
## 5 7.7 19249.23 16.3 40748.37 5.1 12749.490
## 6 7.7 10656.18 15.8 21865.94 2.7 3736.584
## geometry
## 1 MULTIPOLYGON (((7.72106 4.9...
## 2 MULTIPOLYGON (((7.689583 4....
## 3 MULTIPOLYGON (((7.891184 4....
## 4 MULTIPOLYGON (((8.142083 4....
## 5 MULTIPOLYGON (((7.610435 5....
## 6 MULTIPOLYGON (((7.702902 4....View(ssdata) #view data set in a window
ssdata <- st_as_sf(ssdata) #convert object to sf object
class(ssdata)## [1] "sf" "data.frame"str(ssdata) #view structure of an object## Classes 'sf' and 'data.frame': 123 obs. of 36 variables:
## $ NAME_1 : chr "Akwa Ibom" "Akwa Ibom" "Akwa Ibom" "Akwa Ibom" ...
## $ NAME_2 : chr "Abak" "Eastern Obolo" "Eket" "Esit Eket" ...
## $ ID_0 : chr "163" "163" "163" "163" ...
## $ ISO : chr "NGA" "NGA" "NGA" "NGA" ...
## $ NAME_0 : chr "Nigeria" "Nigeria" "Nigeria" "Nigeria" ...
## $ ID_1 : chr "3" "3" "3" "3" ...
## $ ID_2 : chr "39" "40" "41" "42" ...
## $ TYPE_2 : chr "Local Authority" "Local Authority" "Local Authority" "Local Authority" ...
## $ ENGTYPE_2 : chr "Local Authority" "Local Authority" "Local Authority" "Local Authority" ...
## $ NL_NAME_2 : chr NA NA NA NA ...
## $ VARNAME_2 : chr NA NA NA NA ...
## $ OBJECTID : num 39 40 41 42 43 44 45 46 47 48 ...
## $ LGA : chr "Abak" "Eastern Obolo" "Eket" "Esit - Eket" ...
## $ AREA (km2) : num 167 156 210 148 298 ...
## $ POP : num 154945 125525 314900 93225 249990 ...
## $ POP_DEN : num 928 802 1502 629 839 ...
## $ PUBLIC_HF : num 16 12 19 15 21 14 23 20 18 21 ...
## $ PHF_COV : num 10.33 9.56 6.03 16.09 8.4 ...
## $ HIV_PREV : num 5.5 5.5 5.5 5.5 5.5 5.5 5.5 5.5 5.5 5.5 ...
## $ EST_PLHIV : num 8522 6904 17320 5127 13749 ...
## $ TX_CURR : num 5363 5499 17231 6646 9543 ...
## $ NOR_TXCURR : num 3461 4381 5472 7129 3817 ...
## $ ART_HF : num 4 0 17 0 12 0 5 4 0 0 ...
## $ AHF_COV : num 7.46 0 9.87 0 12.57 ...
## $ MPI INDEX : num 0.288 0.308 0.308 0.308 0.288 0.288 0.283 0.308 0.283 0.283 ...
## $ DEPRIVATION: num 71.3 73.2 73.2 73.2 71.3 71.3 69.3 73.2 69.3 69.3 ...
## $ POP_POV : num 110476 91884 230507 68241 178243 ...
## $ HEALTH_DEP%: num 42 39.4 39.4 39.4 42 42 40.5 39.4 40.5 40.5 ...
## $ HEALTH_DEP : num 65077 49457 124071 36731 104996 ...
## $ EDUC_DEP% : num 7.7 6.4 6.4 6.4 7.7 7.7 7.1 6.4 7.1 7.1 ...
## $ EDUC_DEP : num 11931 8034 20154 5966 19249 ...
## $ WORK_DEP% : num 16.4 16.6 20 17.4 16.3 15.8 16.6 16.3 17.4 17.4 ...
## $ WORK_DEP : num 25411 20837 62980 16221 40748 ...
## $ SEC_SHOCK% : num 8 15.2 5.6 5.4 5.1 2.7 15.2 5.1 5.4 5.4 ...
## $ SEC_SHOCK : num 12396 19080 17634 5034 12749 ...
## $ geometry :sfc_MULTIPOLYGON of length 123; first list element: List of 1
## ..$ :List of 1
## .. ..$ : num [1:111, 1:2] 7.72 7.72 7.73 7.73 7.73 ...
## ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "names")= chr [1:35] "NAME_1" "NAME_2" "ID_0" "ISO" ...Data Visualisation of the 123 Local Government Areas in South South Nigeria
Map showing the population of South-South Nigeria
describe(rschdata$POP, fast = TRUE)## rschdata$POP
## n missing distinct Info Mean Gmd .05 .10
## 123 0 123 1 246110 140365 97031 121815
## .25 .50 .75 .90 .95
## 154806 215165 284152 442095 526765
##
## lowest : 48200 79201 83066 91670 93225, highest: 564498 574199 620174 708008 776979options(scipen=100000) #before a plot sets legend numbers)
ggplot(rschdata, aes(POP, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1) + scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Population") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
ggplot(rschdata, aes(POP, fill = factor(NAME_1))) + geom_histogram() + scale_fill_viridis_d(option = "viridis") + labs(x = "Population", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
tm_shape(ssdata) + tm_fill (col = "POP", palette = "-viridis", title = "Population")
Map showing the population density of South-South Nigeria see if u can change legend scale to 500
describe(rschdata$POP_DEN, fast = TRUE)## rschdata$POP_DEN
## n missing distinct Info Mean Gmd .05 .10
## 123 0 123 1 799.8 805.5 81.1 124.7
## .25 .50 .75 .90 .95
## 268.8 489.9 969.5 1547.8 2643.9
##
## lowest : 35.9677 62.8588 73.8839 74.352 74.4682
## highest: 3046.18 3459.24 3904.73 4324.91 5829.43ggplot(rschdata, aes(POP_DEN, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1, binwidth = 150) + scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Population Density") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
ggplot(rschdata, aes(POP_DEN, fill = factor(NAME_1))) + geom_histogram() + scale_fill_viridis_d(option = "viridis") + labs(x = "Population Density", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
tm_shape(ssdata) + tm_fill (col = "POP_DEN", palette = "-viridis", title = "Population Density")
Map showing the distribution of Public Health facilities in South-South
Nigeria
describe(rschdata$PUBLIC_HF, fast = TRUE)## rschdata$PUBLIC_HF
## n missing distinct Info Mean Gmd .05 .10
## 123 0 46 0.998 28.41 17.69 12.0 14.0
## .25 .50 .75 .90 .95
## 17.0 22.0 34.5 51.6 64.3
##
## lowest : 6 8 9 11 12, highest: 65 67 73 92 120ggplot(rschdata, aes(PUBLIC_HF, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1, binwidth = 3.5) + scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Public Health Facility") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
ggplot(rschdata, aes(PUBLIC_HF, fill = factor(NAME_1))) + geom_histogram() + scale_fill_viridis_d(option = "viridis") + labs(x = "Public Health Facility", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
tm_shape(ssdata) + tm_fill (col = "PUBLIC_HF", palette = "-viridis", title = "Public HF")+ tm_layout(legend.position = c(0.45, 0.45))
Map showing the coverage of Public Health facilities per 100,000 persons in South-South Nigeria
describe(rschdata$PHF_COV, fast = TRUE)## rschdata$PHF_COV
## n missing distinct Info Mean Gmd .05 .10
## 123 0 123 1 13.99 10.22 3.833 5.615
## .25 .50 .75 .90 .95
## 7.608 10.356 16.622 28.673 33.709
##
## lowest : 2.27314 2.78856 3.48071 3.57345 3.64434
## highest: 38.31 38.3735 43.9378 48.1705 72.6141ggplot(rschdata, aes(PHF_COV, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1, binwidth = 2.1) + scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Public Health Facility Coverage") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
ggplot(rschdata, aes(PHF_COV, fill = factor(NAME_1))) + geom_histogram() + scale_fill_viridis_d(option = "viridis") + labs(x = "Public Health Facility Coverage", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
tm_shape(ssdata) + tm_fill (col = "PHF_COV", palette = "-viridis", title = "Public HF Coverage")
Map showing the estimated people living with HIV in South-South Nigeria
describe(rschdata$EST_PLHIV, fast = TRUE)## rschdata$EST_PLHIV
## n missing distinct Info Mean Gmd .05 .10
## 123 0 123 1 7578 5389 2347 2909
## .25 .50 .75 .90 .95
## 3717 6057 9634 14703 18687
##
## lowest : 964 1578.25 1833.4 1853.79 2014.36
## highest: 20060.3 21450.9 21819.6 22440.6 26904.3ggplot(rschdata, aes(EST_PLHIV, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1, binwidth = 850) + scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Estimated PLHIV") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
ggplot(rschdata, aes(EST_PLHIV, fill = factor(NAME_1))) + geom_histogram() + scale_fill_viridis_d(option = "viridis") + labs(x = "Estimated PLHIV", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
tm_shape(ssdata) + tm_fill (col = "EST_PLHIV", palette = "-viridis", title = "Estimated PLHIV")
Map showing the distribution of people currently on HIV treatment (TX_CURR) in South-South Nigeria
describe(rschdata$TX_CURR, fast = TRUE)## rschdata$TX_CURR
## n missing distinct Info Mean Gmd .05 .10
## 123 0 111 0.999 4970 5822 0 15
## .25 .50 .75 .90 .95
## 849 2990 6436 11847 15320
##
## lowest : 0 75 96 132 278, highest: 16045 17231 20236 21205 60049ggplot(rschdata, aes(TX_CURR, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1) + scale_fill_viridis_d(option = "viridis") + theme_bw()## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
ggplot(rschdata, aes(TX_CURR, fill = factor(NAME_1))) + geom_histogram() + scale_fill_viridis_d(option = "viridis") + theme_bw()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(ssdata) + geom_sf(aes(fill = TX_CURR)) + scale_fill_viridis(direction = -1) + theme_bw()
Map showing the distribution of people currently on HIV treatment (NORTX_CURR) in South-South Nigeria
describe(rschdata$NOR_TXCURR, fast = TRUE)## rschdata$NOR_TXCURR
## n missing distinct Info Mean Gmd .05 .10
## 123 0 111 0.999 2185 2505 0.000 6.776
## .25 .50 .75 .90 .95
## 399.068 1444.339 3171.242 5384.625 7117.494
##
## lowest : 0 33.8798 82.4284 115.571 125.718
## highest: 8420.32 8481.4 10005.2 10895.6 13298.8ggplot(rschdata, aes(EST_PLHIV, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1, binwidth = 450) + scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Noramlised TX_CURR") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
ggplot(rschdata, aes(EST_PLHIV, fill = factor(NAME_1))) + geom_histogram() + scale_fill_viridis_d(option = "viridis") + labs(x = "Noramlised TX_CURR", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(ssdata) + geom_sf(aes(fill = NOR_TXCURR)) + scale_fill_viridis(direction = -1) + theme_bw()
Map showing the Public Health facilities offering HIV services in South-South Nigeria
describe(rschdata$ART_HF, fast = TRUE)## rschdata$ART_HF
## n missing distinct Info Mean Gmd .05 .10
## 123 0 36 0.997 12.42 11.24 0.0 2.0
## .25 .50 .75 .90 .95
## 5.0 10.0 17.0 27.0 31.8
##
## lowest : 0 1 2 3 4, highest: 35 43 47 51 57ggplot(rschdata, aes(ART_HF, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1, binwidth = 2.2) + scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "ART Health Facility") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
ggplot(rschdata, aes(ART_HF, fill = factor(NAME_1))) + geom_histogram() + scale_fill_viridis_d(option = "viridis") + labs(x = "ART Health Facility", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
tm_shape(ssdata) + tm_fill (col = "ART_HF", palette = "-viridis", title = "ART HF")
Map showing the distribution of Public Health facilities offering HIV services per 10,000 persons on treatment in South-South Nigeria
describe(rschdata$AHF_COV, fast = TRUE)## rschdata$AHF_COV
## n missing distinct Info Mean Gmd .05 .10
## 123 0 105 0.996 97.82 150.2 0.000 0.000
## .25 .50 .75 .90 .95
## 5.154 15.582 74.008 310.402 453.638
##
## lowest : 0 1.48251 2.0024 2.57467 2.66449
## highest: 497.238 511.073 569.106 1066.67 1458.33ggplot(rschdata, aes(AHF_COV, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1, binwidth = 23) + scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "ART Health Facility Coverage") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
ggplot(rschdata, aes(AHF_COV, fill = factor(NAME_1))) + geom_histogram() + scale_fill_viridis_d(option = "viridis") + labs(x = "ART Health Facility Coverage", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
tm_shape(ssdata) + tm_fill (col = "AHF_COV", palette = "-viridis", title = "ART HF Coverage")
Map showing the distribution of population poverty in South-South Nigeria
describe(rschdata$POP_POV, fast = TRUE)## rschdata$POP_POV
## n missing distinct Info Mean Gmd .05 .10
## 123 0 123 1 152156 99149 52314 56797
## .25 .50 .75 .90 .95
## 91804 124811 184881 264225 363202
##
## lowest : 31571 34929.3 36120.7 43609.6 47691.1
## highest: 372700 398536 406214 490204 661986ggplot(rschdata, aes(POP_POV, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1) + scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Population Poverty") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
ggplot(rschdata, aes(POP_POV, fill = factor(NAME_1))) + geom_histogram() + scale_fill_viridis_d(option = "viridis") + labs(x = "Population Poverty", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
tm_shape(ssdata) + tm_fill (col = "POP_POV", palette = "-viridis", title = "Population Poverty")
Map showing the distribution of health deprivation in South-South Nigeria
describe(rschdata$HEALTH_DEP, fast = TRUE)## rschdata$HEALTH_DEP
## n missing distinct Info Mean Gmd .05 .10
## 123 0 123 1 88015 48501 36758 44653
## .25 .50 .75 .90 .95
## 55800 74922 104984 160923 186055
##
## lowest : 18316 28990 30892.8 31205.2 33939.8
## highest: 199268 200970 214446 235666 247803ggplot(rschdata, aes(HEALTH_DEP, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1) + scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Health Deprivation") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
ggplot(rschdata, aes(HEALTH_DEP, fill = factor(NAME_1))) + geom_histogram() + scale_fill_viridis_d(option = "viridis") + labs(x = "Health Deprivation", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
tm_shape(ssdata) + tm_fill (col = "HEALTH_DEP", palette = "-viridis", title = "Health Deprivation")
Map showing the distribution of education deprivation in South-South Nigeria
attach(rschdata)
plot(EDUC_DEP,TX_CURR, type="p")
describe(rschdata$EDUC_DEP, fast = TRUE)## rschdata$EDUC_DEP
## n missing distinct Info Mean Gmd .05 .10
## 123 0 123 1 21360 14958 7052 8071
## .25 .50 .75 .90 .95
## 11487 17207 26599 41255 54597
##
## lowest : 3663.2 5068.86 5430.38 5455.06 5500.2
## highest: 58282.2 60251.8 60645.8 67416.4 79832.2ggplot(rschdata, aes(EDUC_DEP, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1) + scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Education Deprivation") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
ggplot(rschdata, aes(EDUC_DEP, fill = factor(NAME_1))) + geom_histogram() + scale_fill_viridis_d(option = "viridis") + labs(x = "Education Deprivation", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
tm_shape(ssdata) + tm_fill (col = "EDUC_DEP", palette = "-viridis", title = "Education Deprivation")
Map showing the distribution of work deprivation in South-South Nigeria
describe(rschdata$WORK_DEP, fast = TRUE)## rschdata$WORK_DEP
## n missing distinct Info Mean Gmd .05 .10
## 123 0 123 1 38962 25693 12651 15765
## .25 .50 .75 .90 .95
## 21732 32962 47730 68826 87777
##
## lowest : 8531.4 9794.07 9954.53 12125.3 12404.9
## highest: 93594.4 98222.7 98783.4 141602 155396ggplot(rschdata, aes(WORK_DEP, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1) + scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Work Deprivation") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
ggplot(rschdata, aes(WORK_DEP, fill = factor(NAME_1))) + geom_histogram() + scale_fill_viridis_d(option = "viridis") + labs(x = "Work Deprivation", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
tm_shape(ssdata) + tm_fill (col = "WORK_DEP", palette = "-viridis", title = "Work Deprivation")
Map showing the distribution of security shock in South-South Nigeria
describe(rschdata$SEC_SHOCK, fast = TRUE)## rschdata$SEC_SHOCK
## n missing distinct Info Mean Gmd .05 .10
## 123 0 123 1 17555 10894 5666 7366
## .25 .50 .75 .90 .95
## 10583 15130 22254 31170 40245
##
## lowest : 2138.43 3736.58 3826.44 4331.96 4410.4
## highest: 41934.5 43152.6 43510.8 43539.6 49079.7var(rschdata$SEC_SHOCK)## [1] 101242585ggplot(rschdata, aes(SEC_SHOCK, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1) + scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Security Shock") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
ggplot(rschdata, aes(SEC_SHOCK, fill = factor(NAME_1))) + geom_histogram() + scale_fill_viridis_d(option = "viridis") + labs(x = "Security Shock", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
tm_shape(ssdata) + tm_fill (col = "SEC_SHOCK", palette = "viridis", title = "SECURITY_SHOCK") 
Aspatial bivariate analysis using Scatter plot and matrix
Scatter plot
rschvar <- rschdata[c('NOR_TXCURR', 'POP_DEN', 'PHF_COV', 'AHF_COV', 'POP_POV', 'HEALTH_DEP', 'EDUC_DEP', 'WORK_DEP', 'SEC_SHOCK')]
C <- cor(rschvar, method = "pearson", use ="complete.obs")
corrplot.mixed(C, tl.cex = 0.8, tl.col="black")
Cc <- cor(rschvar, method = "spearman", use ="complete.obs")
corrplot.mixed(Cc, tl.cex = 0.8, tl.col="black")
pairs(rschvar)
Scatter Plot Matrix
rcorr(as.matrix(rschvar), type ="pearson")## NOR_TXCURR POP_DEN PHF_COV AHF_COV POP_POV HEALTH_DEP EDUC_DEP
## NOR_TXCURR 1.00 0.24 0.09 -0.33 -0.06 -0.08 -0.26
## POP_DEN 0.24 1.00 -0.31 -0.20 0.26 0.40 0.15
## PHF_COV 0.09 -0.31 1.00 0.08 -0.33 -0.44 -0.31
## AHF_COV -0.33 -0.20 0.08 1.00 -0.05 -0.14 -0.04
## POP_POV -0.06 0.26 -0.33 -0.05 1.00 0.83 0.35
## HEALTH_DEP -0.08 0.40 -0.44 -0.14 0.83 1.00 0.72
## EDUC_DEP -0.26 0.15 -0.31 -0.04 0.35 0.72 1.00
## WORK_DEP 0.00 0.41 -0.42 -0.09 0.77 0.89 0.62
## SEC_SHOCK -0.27 0.20 -0.32 -0.03 0.70 0.78 0.60
## WORK_DEP SEC_SHOCK
## NOR_TXCURR 0.00 -0.27
## POP_DEN 0.41 0.20
## PHF_COV -0.42 -0.32
## AHF_COV -0.09 -0.03
## POP_POV 0.77 0.70
## HEALTH_DEP 0.89 0.78
## EDUC_DEP 0.62 0.60
## WORK_DEP 1.00 0.70
## SEC_SHOCK 0.70 1.00
##
## n= 123
##
##
## P
## NOR_TXCURR POP_DEN PHF_COV AHF_COV POP_POV HEALTH_DEP EDUC_DEP
## NOR_TXCURR 0.0087 0.2981 0.0002 0.5160 0.3743 0.0031
## POP_DEN 0.0087 0.0004 0.0274 0.0034 0.0000 0.0912
## PHF_COV 0.2981 0.0004 0.3636 0.0002 0.0000 0.0004
## AHF_COV 0.0002 0.0274 0.3636 0.5887 0.1249 0.6444
## POP_POV 0.5160 0.0034 0.0002 0.5887 0.0000 0.0000
## HEALTH_DEP 0.3743 0.0000 0.0000 0.1249 0.0000 0.0000
## EDUC_DEP 0.0031 0.0912 0.0004 0.6444 0.0000 0.0000
## WORK_DEP 0.9903 0.0000 0.0000 0.3097 0.0000 0.0000 0.0000
## SEC_SHOCK 0.0023 0.0246 0.0003 0.7643 0.0000 0.0000 0.0000
## WORK_DEP SEC_SHOCK
## NOR_TXCURR 0.9903 0.0023
## POP_DEN 0.0000 0.0246
## PHF_COV 0.0000 0.0003
## AHF_COV 0.3097 0.7643
## POP_POV 0.0000 0.0000
## HEALTH_DEP 0.0000 0.0000
## EDUC_DEP 0.0000 0.0000
## WORK_DEP 0.0000
## SEC_SHOCK 0.0000rcorr(as.matrix(rschvar), type ="spearman")## NOR_TXCURR POP_DEN PHF_COV AHF_COV POP_POV HEALTH_DEP EDUC_DEP
## NOR_TXCURR 1.00 0.44 -0.17 -0.32 0.07 0.00 -0.28
## POP_DEN 0.44 1.00 -0.63 -0.30 0.37 0.39 0.05
## PHF_COV -0.17 -0.63 1.00 0.20 -0.51 -0.62 -0.33
## AHF_COV -0.32 -0.30 0.20 1.00 0.03 0.02 0.15
## POP_POV 0.07 0.37 -0.51 0.03 1.00 0.84 0.31
## HEALTH_DEP 0.00 0.39 -0.62 0.02 0.84 1.00 0.69
## EDUC_DEP -0.28 0.05 -0.33 0.15 0.31 0.69 1.00
## WORK_DEP -0.01 0.41 -0.61 -0.10 0.71 0.86 0.61
## SEC_SHOCK -0.26 0.12 -0.42 0.09 0.62 0.76 0.63
## WORK_DEP SEC_SHOCK
## NOR_TXCURR -0.01 -0.26
## POP_DEN 0.41 0.12
## PHF_COV -0.61 -0.42
## AHF_COV -0.10 0.09
## POP_POV 0.71 0.62
## HEALTH_DEP 0.86 0.76
## EDUC_DEP 0.61 0.63
## WORK_DEP 1.00 0.66
## SEC_SHOCK 0.66 1.00
##
## n= 123
##
##
## P
## NOR_TXCURR POP_DEN PHF_COV AHF_COV POP_POV HEALTH_DEP EDUC_DEP
## NOR_TXCURR 0.0000 0.0537 0.0003 0.4599 0.9980 0.0020
## POP_DEN 0.0000 0.0000 0.0006 0.0000 0.0000 0.5998
## PHF_COV 0.0537 0.0000 0.0303 0.0000 0.0000 0.0002
## AHF_COV 0.0003 0.0006 0.0303 0.7512 0.8386 0.1093
## POP_POV 0.4599 0.0000 0.0000 0.7512 0.0000 0.0005
## HEALTH_DEP 0.9980 0.0000 0.0000 0.8386 0.0000 0.0000
## EDUC_DEP 0.0020 0.5998 0.0002 0.1093 0.0005 0.0000
## WORK_DEP 0.8955 0.0000 0.0000 0.2823 0.0000 0.0000 0.0000
## SEC_SHOCK 0.0038 0.1715 0.0000 0.3035 0.0000 0.0000 0.0000
## WORK_DEP SEC_SHOCK
## NOR_TXCURR 0.8955 0.0038
## POP_DEN 0.0000 0.1715
## PHF_COV 0.0000 0.0000
## AHF_COV 0.2823 0.3035
## POP_POV 0.0000 0.0000
## HEALTH_DEP 0.0000 0.0000
## EDUC_DEP 0.0000 0.0000
## WORK_DEP 0.0000
## SEC_SHOCK 0.0000##Near Neighbour Analysis Centroids
ssdatacast <- st_cast(ssdata, "POLYGON") #change multipolygon to polygon in .shf data set## Warning in st_cast.sf(ssdata, "POLYGON"): repeating attributes for all
## sub-geometries for which they may not be constanthead(ssdatacast)## Simple feature collection with 6 features and 35 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 7.555695 ymin: 4.472917 xmax: 7.828423 ymax: 5.099989
## Geodetic CRS: WGS 84
## NAME_1 NAME_2 ID_0 ISO NAME_0 ID_1 ID_2 TYPE_2
## 1 Akwa Ibom Abak 163 NGA Nigeria 3 39 Local Authority
## 2 Akwa Ibom Eastern Obolo 163 NGA Nigeria 3 40 Local Authority
## 2.1 Akwa Ibom Eastern Obolo 163 NGA Nigeria 3 40 Local Authority
## 2.2 Akwa Ibom Eastern Obolo 163 NGA Nigeria 3 40 Local Authority
## 2.3 Akwa Ibom Eastern Obolo 163 NGA Nigeria 3 40 Local Authority
## 2.4 Akwa Ibom Eastern Obolo 163 NGA Nigeria 3 40 Local Authority
## ENGTYPE_2 NL_NAME_2 VARNAME_2 OBJECTID LGA AREA (km2)
## 1 Local Authority <NA> <NA> 39 Abak 167.0
## 2 Local Authority <NA> <NA> 40 Eastern Obolo 156.5
## 2.1 Local Authority <NA> <NA> 40 Eastern Obolo 156.5
## 2.2 Local Authority <NA> <NA> 40 Eastern Obolo 156.5
## 2.3 Local Authority <NA> <NA> 40 Eastern Obolo 156.5
## 2.4 Local Authority <NA> <NA> 40 Eastern Obolo 156.5
## POP POP_DEN PUBLIC_HF PHF_COV HIV_PREV EST_PLHIV TX_CURR NOR_TXCURR
## 1 154945 927.8144 16 10.326245 5.5 8521.975 5363 3461.228
## 2 125525 802.0767 12 9.559849 5.5 6903.875 5499 4380.801
## 2.1 125525 802.0767 12 9.559849 5.5 6903.875 5499 4380.801
## 2.2 125525 802.0767 12 9.559849 5.5 6903.875 5499 4380.801
## 2.3 125525 802.0767 12 9.559849 5.5 6903.875 5499 4380.801
## 2.4 125525 802.0767 12 9.559849 5.5 6903.875 5499 4380.801
## ART_HF AHF_COV MPI INDEX DEPRIVATION POP_POV HEALTH_DEP% HEALTH_DEP
## 1 4 7.458512 0.288 71.3 110475.8 42.0 65076.90
## 2 0 0.000000 0.308 73.2 91884.3 39.4 49456.85
## 2.1 0 0.000000 0.308 73.2 91884.3 39.4 49456.85
## 2.2 0 0.000000 0.308 73.2 91884.3 39.4 49456.85
## 2.3 0 0.000000 0.308 73.2 91884.3 39.4 49456.85
## 2.4 0 0.000000 0.308 73.2 91884.3 39.4 49456.85
## EDUC_DEP% EDUC_DEP WORK_DEP% WORK_DEP SEC_SHOCK% SEC_SHOCK
## 1 7.7 11930.77 16.4 25410.98 8.0 12395.6
## 2 6.4 8033.60 16.6 20837.15 15.2 19079.8
## 2.1 6.4 8033.60 16.6 20837.15 15.2 19079.8
## 2.2 6.4 8033.60 16.6 20837.15 15.2 19079.8
## 2.3 6.4 8033.60 16.6 20837.15 15.2 19079.8
## 2.4 6.4 8033.60 16.6 20837.15 15.2 19079.8
## geometry
## 1 POLYGON ((7.72106 4.965561,...
## 2 POLYGON ((7.689583 4.516528...
## 2.1 POLYGON ((7.619028 4.474583...
## 2.2 POLYGON ((7.645416 4.525392...
## 2.3 POLYGON ((7.67736 4.519306,...
## 2.4 POLYGON ((7.555695 4.504077...st_centroid(ssdata)## Warning: st_centroid assumes attributes are constant over geometries## Simple feature collection with 123 features and 35 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 5.236942 ymin: 4.462525 xmax: 9.28852 ymax: 7.342454
## Geodetic CRS: WGS 84
## First 10 features:
## NAME_1 NAME_2 ID_0 ISO NAME_0 ID_1 ID_2 TYPE_2
## 1 Akwa Ibom Abak 163 NGA Nigeria 3 39 Local Authority
## 2 Akwa Ibom Eastern Obolo 163 NGA Nigeria 3 40 Local Authority
## 3 Akwa Ibom Eket 163 NGA Nigeria 3 41 Local Authority
## 4 Akwa Ibom Esit Eket 163 NGA Nigeria 3 42 Local Authority
## 5 Akwa Ibom Essien-U 163 NGA Nigeria 3 43 Local Authority
## 6 Akwa Ibom EtimEkpo 163 NGA Nigeria 3 44 Local Authority
## 7 Akwa Ibom Etinan 163 NGA Nigeria 3 45 Local Authority
## 8 Akwa Ibom Ibeno 163 NGA Nigeria 3 46 Local Authority
## 9 Akwa Ibom Ibesikpo Asutan 163 NGA Nigeria 3 47 Local Authority
## 10 Akwa Ibom Ibiono Ibom 163 NGA Nigeria 3 48 Local Authority
## ENGTYPE_2 NL_NAME_2 VARNAME_2 OBJECTID LGA AREA (km2)
## 1 Local Authority <NA> <NA> 39 Abak 167.0
## 2 Local Authority <NA> <NA> 40 Eastern Obolo 156.5
## 3 Local Authority <NA> <NA> 41 Eket 209.7
## 4 Local Authority <NA> <NA> 42 Esit - Eket 148.2
## 5 Local Authority <NA> Essien Udim 43 Essien Udim 298.1
## 6 Local Authority <NA> Etim Ekpo 44 Etim Ekpo 200.4
## 7 Local Authority <NA> <NA> 45 Etinan 157.2
## 8 Local Authority <NA> <NA> 46 Ibeno 271.2
## 9 Local Authority <NA> <NA> 47 Ibesikpo Asutan 175.5
## 10 Local Authority <NA> <NA> 48 Ibiono Ibom 336.0
## POP POP_DEN PUBLIC_HF PHF_COV HIV_PREV EST_PLHIV TX_CURR NOR_TXCURR
## 1 154945 927.8144 16 10.326245 5.5 8521.975 5363 3461.228
## 2 125525 802.0767 12 9.559849 5.5 6903.875 5499 4380.801
## 3 314900 1501.6691 19 6.033661 5.5 17319.500 17231 5471.896
## 4 93225 629.0486 15 16.090105 5.5 5127.375 6646 7128.989
## 5 249990 838.6112 21 8.400336 5.5 13749.450 9543 3817.353
## 6 138392 690.5788 14 10.116192 5.5 7611.560 4469 3229.233
## 7 275885 1754.9936 23 8.336807 5.5 15173.675 5376 1948.638
## 8 110123 406.0583 20 18.161510 5.5 6056.765 11018 10005.176
## 9 217014 1236.5470 18 8.294396 5.5 11935.770 6594 3038.514
## 10 279438 831.6607 21 7.515084 5.5 15369.090 8008 2865.752
## ART_HF AHF_COV MPI INDEX DEPRIVATION POP_POV HEALTH_DEP% HEALTH_DEP
## 1 4 7.458512 0.288 71.3 110475.79 42.0 65076.90
## 2 0 0.000000 0.308 73.2 91884.30 39.4 49456.85
## 3 17 9.865939 0.308 73.2 230506.80 39.4 124070.60
## 4 0 0.000000 0.308 73.2 68240.70 39.4 36730.65
## 5 12 12.574662 0.288 71.3 178242.87 42.0 104995.80
## 6 0 0.000000 0.288 71.3 98673.50 42.0 58124.64
## 7 5 9.300595 0.283 69.3 191188.30 40.5 111733.43
## 8 4 3.630423 0.308 73.2 80610.04 39.4 43388.46
## 9 0 0.000000 0.283 69.3 150390.70 40.5 87890.67
## 10 0 0.000000 0.283 69.3 193650.53 40.5 113172.39
## EDUC_DEP% EDUC_DEP WORK_DEP% WORK_DEP SEC_SHOCK% SEC_SHOCK
## 1 7.7 11930.765 16.4 25410.98 8.0 12395.600
## 2 6.4 8033.600 16.6 20837.15 15.2 19079.800
## 3 6.4 20153.600 20.0 62980.00 5.6 17634.400
## 4 6.4 5966.400 17.4 16221.15 5.4 5034.150
## 5 7.7 19249.230 16.3 40748.37 5.1 12749.490
## 6 7.7 10656.184 15.8 21865.94 2.7 3736.584
## 7 7.1 19587.835 16.6 45796.91 15.2 41934.520
## 8 6.4 7047.872 16.3 17950.05 5.1 5616.273
## 9 7.1 15407.994 17.4 37760.44 5.4 11718.756
## 10 7.1 19840.098 17.4 48622.21 5.4 15089.652
## geometry
## 1 POINT (7.760623 5.006096)
## 2 POINT (7.696764 4.516807)
## 3 POINT (7.94293 4.628248)
## 4 POINT (8.067648 4.633393)
## 5 POINT (7.638981 5.099096)
## 6 POINT (7.590369 4.960887)
## 7 POINT (7.834807 4.801225)
## 8 POINT (8.09915 4.561915)
## 9 POINT (7.940853 4.911513)
## 10 POINT (7.881161 5.199924)plot(st_geometry(ssdata))
plot(st_geometry(st_centroid(ssdata)), pch = 3, cex = 0.5, col = 'red', add=TRUE)## Warning: st_centroid assumes attributes are constant over geometries
Near Neighbour Analysis
ID_2 = row.names(ssdatacast)
ssdatacast$ID_2 = ID_2
sscastnb1 <- poly2nb(as(ssdata, "Spatial")) # creating adjacency matrix (neighborhood graph)
summary(sscastnb1)## Neighbour list object:
## Number of regions: 123
## Number of nonzero links: 630
## Percentage nonzero weights: 4.164188
## Average number of links: 5.121951
## Link number distribution:
##
## 1 2 3 4 5 6 7 8 9 10 11
## 2 5 12 34 24 18 17 5 2 3 1
## 2 least connected regions:
## 43 119 with 1 link
## 1 most connected region:
## 95 with 11 linkshead(sscastnb1)## [[1]]
## [1] 5 6 12 18 26 28 31
##
## [[2]]
## [1] 8 13 18 24 120
##
## [[3]]
## [1] 4 7 8 21 24
##
## [[4]]
## [1] 3 8 17 21 23 30
##
## [[5]]
## [1] 1 6 11 12 14 22
##
## [[6]]
## [1] 1 5 11 28plot(st_geometry(ssdata), border="grey60", axes = TRUE)
plot(sscastnb1, st_coordinates(st_centroid(ssdata)), pch = 19, cex = 0.01, col="red", add=TRUE)## Warning: st_centroid assumes attributes are constant over geometries
head(st_coordinates(st_centroid(ssdata)))## Warning: st_centroid assumes attributes are constant over geometries## X Y
## [1,] 7.760623 5.006096
## [2,] 7.696764 4.516807
## [3,] 7.942930 4.628248
## [4,] 8.067648 4.633393
## [5,] 7.638981 5.099096
## [6,] 7.590369 4.960887Near Neighbor 2
sscastcoords <- st_coordinates(st_centroid(st_geometry(ssdata)))
sscastnbsf = nb2lines(sscastnb1, coords=sscastcoords, proj4string=CRS("+proj=longlat +datum=WGS84"), as_sf=T) # weights are binary
summary(sscastnbsf)## i j i_ID j_ID
## Min. : 1.00 Min. : 1.00 Length:630 Length:630
## 1st Qu.: 30.25 1st Qu.: 30.25 Class :character Class :character
## Median : 63.00 Median : 63.00 Mode :character Mode :character
## Mean : 61.89 Mean : 61.89
## 3rd Qu.: 93.00 3rd Qu.: 93.00
## Max. :123.00 Max. :123.00
## wt geometry
## Min. :1 LINESTRING :630
## 1st Qu.:1 epsg:NA : 0
## Median :1 +proj=long...: 0
## Mean :1
## 3rd Qu.:1
## Max. :1ggplot(ssdata) + geom_sf(fill = 'NA', color = 'grey', border='grey60') + geom_sf(data = sscastnbsf, pch = 3, cex = 0.5, col="red") + theme_minimal() + ylab("Latitude") + xlab("Longitude")## Warning in layer_sf(geom = GeomSf, data = data, mapping = mapping, stat = stat,
## : Ignoring unknown parameters: `border`
Moran I Test
moran(ssdata$NOR_TXCURR, nb2listw(sscastnb1, style="W", zero.policy=TRUE), length(sscastnb1), Szero(nb2listw(sscastnb1, style="W", zero.policy=TRUE)), zero.policy = TRUE)## $I
## [1] 0.3333579
##
## $K
## [1] 6.833905moran.test(ssdata$NOR_TXCURR, nb2listw(sscastnb1, style="W", zero.policy=TRUE), na.action=na.pass, zero.policy = TRUE, alternative = "two.sided") #zero.policy for missing values, na.pass puts a 0 instead of deleting observation)##
## Moran I test under randomisation
##
## data: ssdata$NOR_TXCURR
## weights: nb2listw(sscastnb1, style = "W", zero.policy = TRUE)
##
## Moran I statistic standard deviate = 6.024, p-value = 0.000000001701
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.333357854 -0.008196721 0.003214733moran.mc(ssdata$NOR_TXCURR, nb2listw(sscastnb1, style = "W", zero.policy = TRUE), zero.policy = TRUE, nsim = 999, alternative = "two.sided")##
## Monte-Carlo simulation of Moran I
##
## data: ssdata$NOR_TXCURR
## weights: nb2listw(sscastnb1, style = "W", zero.policy = TRUE)
## number of simulations + 1: 1000
##
## statistic = 0.33336, observed rank = 1000, p-value <
## 0.00000000000000022
## alternative hypothesis: two.sidedmoran.plot(ssdata$NOR_TXCURR, nb2listw(sscastnb1, zero.policy=TRUE))
moran.plot(as.vector(scale(ssdata$NOR_TXCURR)), nb2listw(sscastnb1, zero.policy=TRUE), xlim=c(-0.8, 8), ylim=c(-0.8,2.8), pch=1, ylab = 'Spatially Lagged NOR_TXCURR', xlab = 'Normalised TX_CURR')
sslmoran <-localmoran(ssdata$NOR_TXCURR, nb2listw(sscastnb1, style="W", zero.policy=TRUE), zero.policy = TRUE, na.action=na.pass)
print(sslmoran)## Ii E.Ii Var.Ii Z.Ii Pr(z != E(Ii))
## 1 0.262799440 -0.002141948541 0.03569417432 1.4023330 0.1608158453402
## 2 0.453875507 -0.006340871331 0.14987250485 1.1887789 0.2345266842014
## 3 1.606879977 -0.014208229516 0.33316605836 2.8085119 0.0049771040414
## 4 3.551788973 -0.032145961061 0.61145255071 4.5833048 0.0000045768422
## 5 0.339062502 -0.003504176532 0.06862588206 1.3076790 0.1909822112521
## 6 0.153755786 -0.001433974465 0.04293979441 0.7489160 0.4539078643927
## 7 -0.041213007 -0.000073497388 0.00144433011 -1.0824950 0.2790326542547
## 8 4.543410962 -0.080427934843 1.45351382586 3.8352450 0.0001254392000
## 9 0.129015453 -0.000957989160 0.02276564423 0.8614194 0.3890070965404
## 10 0.088040343 -0.000609407414 0.01448700431 0.7365262 0.4614105042275
## 11 0.110555948 -0.000348649238 0.02125730822 0.7606685 0.4468551106494
## 12 0.156166229 -0.000444631917 0.00742210394 1.8178513 0.0690868588962
## 13 -0.366147633 -0.004179024541 0.09899010613 -1.1504692 0.2499506550092
## 14 0.448276411 -0.030668534366 1.19870075876 0.4374518 0.6617837272740
## 15 -0.030969898 -0.001322395264 0.03960302440 -0.1489786 0.8815705047423
## 16 0.289536895 -0.001975592699 0.05912624955 1.1988559 0.2305839685432
## 17 6.280010348 -0.099786205220 2.13674043913 4.3644643 0.0000127434661
## 18 -0.051854263 -0.000028272201 0.00040952644 -2.5609826 0.0104376569128
## 19 0.034514414 -0.000008231447 0.00024683933 2.1973381 0.0279963073197
## 20 -0.283654201 -0.005814270396 0.17334228101 -0.6673327 0.5045596417553
## 21 -0.008932743 -0.000001068861 0.00001548301 -2.2698919 0.0232141412705
## 22 -0.154286451 -0.000300587030 0.00901117517 -1.6221467 0.1047719322864
## 23 1.861093252 -0.048346784180 0.58722155218 2.4917532 0.0127114316274
## 24 0.943028028 -0.009841586946 0.19151213394 2.1773855 0.0294518151042
## 25 2.931836118 -0.027463530972 0.63532781279 3.7127013 0.0002050588576
## 26 0.461585200 -0.010686159535 0.20776970272 1.0360972 0.3001567987920
## 27 3.364479804 -0.022328773206 0.65463535048 4.1859200 0.0000284013372
## 28 0.087935644 -0.000080568057 0.00241584829 1.7907207 0.0733381308427
## 29 1.006944880 -0.020916896615 0.40247905588 1.6201793 0.1051937713384
## 30 0.923285609 -0.001349087713 0.03204715556 5.1650608 0.0000002403607
## 31 0.399350353 -0.010124858943 0.11410565364 1.2121993 0.2254360931278
## 32 0.419899345 -0.002618498178 0.06212253706 1.6951981 0.0900378429644
## 33 0.511401589 -0.005565571029 0.16596925275 1.2689635 0.2044540765636
## 34 0.465561924 -0.006279078720 0.25159720307 0.9406819 0.3468679126656
## 35 0.455194849 -0.003935072872 0.15804686124 1.1548954 0.2481332545047
## 36 0.362016241 -0.003076381111 0.06027379011 1.4870953 0.1369896382763
## 37 0.290721050 -0.004538262499 0.06544082860 1.1541951 0.2484202099147
## 38 0.495395881 -0.005155176712 0.08564823995 1.7103664 0.0871981361078
## 39 0.031060297 -0.000028177181 0.00055374750 1.3211238 0.1864600830176
## 40 0.222690004 -0.004025627438 0.24454159173 0.4584639 0.6466192282599
## 41 0.064132501 -0.000182334839 0.00304445813 1.1656178 0.2437690159198
## 42 -0.856572069 -0.003895904060 0.09230996323 -2.8064663 0.0050088148951
## 43 -3.069109149 -0.162441410889 16.73466646684 -0.7105373 0.4773710307771
## 44 0.131445909 -0.006279078720 0.25159720307 0.2745743 0.7836433366700
## 45 0.248690770 -0.004169508458 0.12451223711 0.7165964 0.4736231956126
## 46 0.117754167 -0.003280919562 0.07778645034 0.4339696 0.6643105164214
## 47 -1.063438040 -0.051131672413 1.95632676267 -0.7237545 0.4692164739028
## 48 0.037642169 -0.000009794150 0.00023296898 2.4668263 0.0136316451239
## 49 0.113495716 -0.006279078720 0.25159720307 0.2387880 0.8112699620346
## 50 -0.156529041 -0.001264667477 0.02482294154 -0.9854741 0.3243912701523
## 51 0.288744767 -0.003164862206 0.19241952395 0.6654629 0.5057544601694
## 52 0.043856871 -0.002758143368 0.08248198284 0.1623104 0.8710614460200
## 53 0.123082893 -0.000721542122 0.02162170668 0.8419595 0.3998106071470
## 54 -0.384187684 -0.003176207340 0.03604665784 -2.0068068 0.0447702424843
## 55 -0.191306225 -0.001809947673 0.04297489807 -0.9140990 0.3606648056928
## 56 -0.041627846 -0.000049474610 0.00148355158 -1.0794838 0.2803720900510
## 57 0.092412711 -0.003189152856 0.09533005594 0.3096360 0.7568377846068
## 58 0.572555629 -0.004641241212 0.13853373374 1.5507657 0.1209578337871
## 59 0.222094682 -0.004180346625 0.06952044265 0.8581839 0.3907909232812
## 60 0.056295509 -0.000040286231 0.00095824098 1.8198980 0.0687745340556
## 61 0.290816374 -0.006279078720 0.18711220564 0.6868232 0.4921941593385
## 62 0.121474067 -0.000594351236 0.01412929747 1.0269352 0.3044509877977
## 63 0.209191113 -0.000832308087 0.02493815123 1.3299509 0.1835344562680
## 64 0.433641344 -0.004303022046 0.08420293704 1.5092296 0.1312401097096
## 65 -0.183487935 -0.000588773016 0.01764549633 -1.3768758 0.1685506406435
## 66 0.300408454 -0.003508015150 0.10482793414 0.9386760 0.3478971116224
## 67 0.296058013 -0.002451223516 0.07332613800 1.1023725 0.2702997645785
## 68 0.074923887 -0.000339590729 0.00491748363 1.0732795 0.2831457464410
## 69 0.271924710 -0.003025616239 0.05037526962 1.2250267 0.2205651680142
## 70 0.334743568 -0.005204003728 0.08645520982 1.1561558 0.2476174385511
## 71 0.325323129 -0.004986583678 0.14879002095 0.8563168 0.3918225963041
## 72 -0.352626141 -0.003182789262 0.12792895718 -0.9769950 0.3285716145660
## 73 0.384928186 -0.005632420217 0.11006987794 1.1772107 0.2391114269497
## 74 -0.197314166 -0.001062137328 0.03181712275 -1.1002306 0.2712316804429
## 75 -0.250683659 -0.002109326935 0.06312023673 -0.9894001 0.3224674024121
## 76 0.311845498 -0.001879438618 0.02135742610 2.1467160 0.0318159016602
## 77 0.223297680 -0.004521422177 0.07516687934 0.8309538 0.4059997323292
## 78 0.345017813 -0.003549084637 0.10605081830 1.0703582 0.2844581187779
## 79 0.318080212 -0.003640835827 0.08628843937 1.0952245 0.2734182764771
## 80 0.105654545 -0.000517322521 0.00863487784 1.1425662 0.2532187625257
## 81 0.078696762 -0.006279078720 0.12262720822 0.2426622 0.8082671268726
## 82 -0.137694467 -0.000776302934 0.01845140965 -1.0079669 0.3134703310049
## 83 0.565098123 -0.005135797196 0.20602366942 1.2563033 0.2090060386896
## 84 0.089015647 -0.000177007017 0.00530707665 1.2243382 0.2208246627724
## 85 0.287257549 -0.001944426205 0.05819530477 1.1988291 0.2305943995008
## 86 0.142726903 -0.000333170293 0.00792240702 1.6072738 0.1079943173271
## 87 0.578669726 -0.005577320481 0.10899914742 1.7696402 0.0767870976051
## 88 0.564855511 -0.005312316044 0.10384775806 1.7693109 0.0768420109934
## 89 0.520546970 -0.006279078720 0.18711220564 1.2179128 0.2232571261932
## 90 0.625521264 -0.006279078720 0.25159720307 1.2595835 0.2078196617412
## 91 0.328960688 -0.001866956745 0.03112018934 1.8753429 0.0607455715388
## 92 0.576489111 -0.006279078720 0.10420292324 1.8053270 0.0710235207650
## 93 0.363532456 -0.006279078720 0.10420292324 1.1456198 0.2519525106859
## 94 -0.128226605 -0.000336903287 0.01358014341 -1.0974466 0.2724462122745
## 95 0.390899548 -0.003928033488 0.04013432229 1.9708316 0.0487431412103
## 96 0.421166472 -0.006279078720 0.10420292324 1.3241612 0.1854495224579
## 97 0.403632792 -0.004209567905 0.25566805064 0.8065924 0.4199013806112
## 98 0.548449695 -0.005105481269 0.15231949193 1.4183492 0.1560888325270
## 99 0.590218568 -0.004083741089 0.12196150629 1.7017511 0.0888020455516
## 100 0.622149899 -0.006279078720 0.10420292324 1.9467772 0.0515614513102
## 101 0.388036777 -0.004664206378 0.07752948910 1.4103553 0.1584348038729
## 102 0.002580034 -0.000000453478 0.00001359871 0.6997660 0.4840734551855
## 103 0.123357333 -0.003805859225 0.09018458304 0.4234433 0.6719718542450
## 104 0.491844751 -0.006085869673 0.24390289803 1.0082317 0.3133432364167
## 105 -0.210218418 -0.006279078720 0.18711220564 -0.4714655 0.6373083277561
## 106 0.358491796 -0.005535466117 0.33574881421 0.6282416 0.5298457064391
## 107 -0.165031090 -0.002566726942 0.05031413386 -0.7242910 0.4688870768594
## 108 0.020238390 -0.001265094420 0.01830238870 0.1589480 0.8737098788773
## 109 -0.270436816 -0.000990209999 0.01944122076 -1.9324620 0.0533025136999
## 110 -0.071160894 -0.001885458427 0.03142801018 -0.3907696 0.6959675118241
## 111 -0.138727146 -0.000397493033 0.01191512533 -1.2672607 0.2050620754513
## 112 0.229389483 -0.006922429959 0.16352243355 0.5843822 0.5589632305307
## 113 -0.021213163 -0.000009997125 0.00040310319 -1.0560697 0.2909363503425
## 114 0.086941147 -0.017871166983 0.22401497393 0.2214490 0.8247428182119
## 115 -0.348641590 -0.052138394127 0.82531953386 -0.3263761 0.7441398363457
## 116 -0.191330838 -0.002284054991 0.04478576003 -0.8933041 0.3716943550869
## 117 0.060382436 -0.000359717976 0.00855345328 0.6567793 0.5113228454889
## 118 -0.084801368 -0.001533515921 0.03642146236 -0.4363135 0.6626092724731
## 119 0.080282714 -0.001089418899 0.13385254408 0.2224141 0.8239915098185
## 120 -0.167793089 -0.006279078720 0.18711220564 -0.3733870 0.7088604661011
## 121 0.351800804 -0.001293655837 0.03073209448 2.0141640 0.0439923206612
## 122 0.061032925 -0.000099646386 0.00298785855 1.1183884 0.2634011543756
## 123 0.529263270 -0.008144690446 0.19215801501 1.2259559 0.2202152756321
## attr(,"call")
## localmoran(x = ssdata$NOR_TXCURR, listw = nb2listw(sscastnb1,
## style = "W", zero.policy = TRUE), zero.policy = TRUE, na.action = na.pass)
## attr(,"class")
## [1] "localmoran" "matrix" "array"
## attr(,"quadr")
## mean median pysal
## 1 High-High High-High High-High
## 2 High-High High-High High-High
## 3 High-High High-High High-High
## 4 High-High High-High High-High
## 5 High-High High-High High-High
## 6 High-High High-High High-High
## 7 Low-High High-High Low-High
## 8 High-High High-High High-High
## 9 High-High High-High High-High
## 10 High-High High-High High-High
## 11 High-High High-High High-High
## 12 High-High High-High High-High
## 13 Low-High Low-High Low-High
## 14 High-High High-High High-High
## 15 High-Low High-High High-Low
## 16 High-High High-High High-High
## 17 High-High High-High High-High
## 18 Low-High High-High Low-High
## 19 High-High High-High High-High
## 20 Low-High Low-High Low-High
## 21 Low-High High-High Low-High
## 22 Low-High High-High Low-High
## 23 High-High High-High High-High
## 24 High-High High-High High-High
## 25 High-High High-High High-High
## 26 High-High High-High High-High
## 27 High-High High-High High-High
## 28 High-High High-High High-High
## 29 High-High High-High High-High
## 30 High-High High-High High-High
## 31 High-High High-High High-High
## 32 Low-Low Low-Low Low-Low
## 33 Low-Low Low-Low Low-Low
## 34 Low-Low Low-Low Low-Low
## 35 Low-Low Low-Low Low-Low
## 36 Low-Low Low-Low Low-Low
## 37 Low-Low Low-Low Low-Low
## 38 Low-Low Low-Low Low-Low
## 39 Low-Low High-Low Low-Low
## 40 Low-Low Low-High Low-Low
## 41 Low-Low High-Low Low-Low
## 42 Low-High Low-High Low-High
## 43 High-Low High-Low High-Low
## 44 Low-Low Low-High Low-Low
## 45 Low-Low Low-Low Low-Low
## 46 Low-Low Low-High Low-Low
## 47 High-Low High-Low High-Low
## 48 High-High High-High High-High
## 49 Low-Low Low-High Low-Low
## 50 High-Low High-Low High-Low
## 51 Low-Low Low-Low Low-Low
## 52 Low-Low Low-High Low-Low
## 53 Low-Low Low-Low Low-Low
## 54 Low-High Low-High Low-High
## 55 High-Low High-Low High-Low
## 56 High-Low High-Low High-Low
## 57 Low-Low Low-High Low-Low
## 58 Low-Low Low-Low Low-Low
## 59 Low-Low Low-High Low-Low
## 60 Low-Low High-Low Low-Low
## 61 Low-Low Low-Low Low-Low
## 62 Low-Low High-Low Low-Low
## 63 Low-Low Low-Low Low-Low
## 64 Low-Low Low-Low Low-Low
## 65 High-Low High-Low High-Low
## 66 Low-Low Low-Low Low-Low
## 67 Low-Low Low-Low Low-Low
## 68 Low-Low High-Low Low-Low
## 69 Low-Low Low-Low Low-Low
## 70 Low-Low Low-Low Low-Low
## 71 Low-Low Low-Low Low-Low
## 72 High-Low High-Low High-Low
## 73 Low-Low Low-Low Low-Low
## 74 High-Low High-Low High-Low
## 75 High-Low High-Low High-Low
## 76 Low-Low Low-Low Low-Low
## 77 Low-Low Low-High Low-Low
## 78 Low-Low Low-Low Low-Low
## 79 Low-Low Low-Low Low-Low
## 80 Low-Low High-Low Low-Low
## 81 Low-Low Low-High Low-Low
## 82 High-Low High-Low High-Low
## 83 Low-Low Low-Low Low-Low
## 84 Low-Low High-Low Low-Low
## 85 Low-Low Low-Low Low-Low
## 86 Low-Low High-Low Low-Low
## 87 Low-Low Low-Low Low-Low
## 88 Low-Low Low-Low Low-Low
## 89 Low-Low Low-Low Low-Low
## 90 Low-Low Low-Low Low-Low
## 91 Low-Low Low-Low Low-Low
## 92 Low-Low Low-Low Low-Low
## 93 Low-Low Low-Low Low-Low
## 94 High-Low High-Low High-Low
## 95 Low-Low Low-Low Low-Low
## 96 Low-Low Low-Low Low-Low
## 97 Low-Low Low-Low Low-Low
## 98 Low-Low Low-Low Low-Low
## 99 Low-Low Low-Low Low-Low
## 100 Low-Low Low-Low Low-Low
## 101 Low-Low Low-Low Low-Low
## 102 Low-Low High-Low Low-Low
## 103 Low-Low Low-High Low-Low
## 104 Low-Low Low-Low Low-Low
## 105 Low-High Low-High Low-High
## 106 Low-Low Low-Low Low-Low
## 107 Low-High Low-High Low-High
## 108 Low-Low Low-High Low-Low
## 109 Low-High Low-High Low-High
## 110 Low-High Low-High Low-High
## 111 Low-High High-High Low-High
## 112 High-High High-High High-High
## 113 Low-High High-High Low-High
## 114 High-High High-High High-High
## 115 High-Low High-High High-Low
## 116 High-Low High-Low High-Low
## 117 High-High High-High High-High
## 118 High-Low High-High High-Low
## 119 Low-Low Low-High Low-Low
## 120 Low-High Low-High Low-High
## 121 High-High High-High High-High
## 122 High-High High-High High-High
## 123 High-High High-High High-Highhist(sslmoran)
lmoranmap <- cbind(ssdata, sslmoran)
mdf <- cbind(ssdata, sslmoran, (attributes(sslmoran)$quadr))
tm_shape(lmoranmap) + tm_fill(col = "Ii", style ="quantile", midpoint = 0, palette = "viridis", title = "Local Moran Observed Statistic") 
tm_shape(lmoranmap) + tm_fill(col = "E.Ii", style ="quantile", midpoint = 0, palette = "-viridis", title = "Local Moran Expected Statistic") 
tm_shape(lmoranmap) + tm_fill(col = "Z.Ii", style ="quantile", midpoint = 0, palette = "-viridis", title = "Local Moran Z-Score Statistic")
tm_shape(lmoranmap) + tm_fill(col = "Pr.z....E.Ii..", style ="quantile", midpoint = 0, palette = "viridis", title = "Local Moran P-Value")
sscastlistw <- nb2listw(sscastnb1, style="W", zero.policy=TRUE)
geary.test(ssdata$NOR_TXCURR, sscastlistw, zero.policy = TRUE)##
## Geary C test under randomisation
##
## data: ssdata$NOR_TXCURR
## weights: sscastlistw
##
## Geary C statistic standard deviate = 5.3765, p-value = 0.00000003798
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.610042600 1.000000000 0.005260618geary.mc(ssdata$NOR_TXCURR,sscastlistw, nsim=999, zero.policy = TRUE)##
## Monte-Carlo simulation of Geary C
##
## data: ssdata$NOR_TXCURR
## weights: sscastlistw
## number of simulations + 1: 1000
##
## statistic = 0.61004, observed rank = 1, p-value = 0.001
## alternative hypothesis: greaterHotspot analysis
sscastlistb <- nb2listw(sscastnb1, style="B", zero.policy=TRUE)
GGis <- globalG.test(ssdata$NOR_TXCURR, sscastlistb, zero.policy = TRUE, alternative = "two.sided")
print(GGis)##
## Getis-Ord global G statistic
##
## data: ssdata$NOR_TXCURR
## weights: sscastlistb
##
## standard deviate = 4.9519, p-value = 0.0000007348
## alternative hypothesis: two.sided
## sample estimates:
## Global G statistic Expectation Variance
## 0.06349260630 0.04198320672 0.00001886716Gis <- localG(ssdata$NOR_TXCURR, sscastlistw, zero.policy = TRUE)
print(Gis)## [1] 1.4023330 1.1887789 2.8085119 4.5833048 1.3076790 0.7489160
## [7] 1.0824950 3.8352450 0.8614194 0.7365262 0.7606685 1.8178513
## [13] 1.1504692 0.4374518 -0.1489786 1.1988559 4.3644643 2.5609826
## [19] 2.1973381 0.6673327 2.2698919 1.6221467 2.4917532 2.1773855
## [25] 3.7127013 1.0360972 4.1859200 1.7907207 1.6201793 5.1650608
## [31] 1.2121993 -1.6951981 -1.2689635 -0.9406819 -1.1548954 -1.4870953
## [37] -1.1541951 -1.7103664 -1.3211238 -0.4584639 -1.1656178 2.8064663
## [43] -0.7105373 -0.2745743 -0.7165964 -0.4339696 -0.7237545 2.4668263
## [49] -0.2387880 -0.9854741 -0.6654629 -0.1623104 -0.8419595 2.0068068
## [55] -0.9140990 -1.0794838 -0.3096360 -1.5507657 -0.8581839 -1.8198980
## [61] -0.6868232 -1.0269352 -1.3299509 -1.5092296 -1.3768758 -0.9386760
## [67] -1.1023725 -1.0732795 -1.2250267 -1.1561558 -0.8563168 -0.9769950
## [73] -1.1772107 -1.1002306 -0.9894001 -2.1467160 -0.8309538 -1.0703582
## [79] -1.0952245 -1.1425662 -0.2426622 -1.0079669 -1.2563033 -1.2243382
## [85] -1.1988291 -1.6072738 -1.7696402 -1.7693109 -1.2179128 -1.2595835
## [91] -1.8753429 -1.8053270 -1.1456198 -1.0974466 -1.9708316 -1.3241612
## [97] -0.8065924 -1.4183492 -1.7017511 -1.9467772 -1.4103553 -0.6997660
## [103] -0.4234433 -1.0082317 0.4714655 -0.6282416 0.7242910 -0.1589480
## [109] 1.9324620 0.3907696 1.2672607 0.5843822 1.0560697 0.2214490
## [115] -0.3263761 -0.8933041 0.6567793 -0.4363135 -0.2224141 0.3733870
## [121] 2.0141640 1.1183884 1.2259559
## attr(,"internals")
## Gi E(Gi) V(Gi) Z(Gi) Pr(z != E(Gi))
## [1,] 0.0130738227 0.008196721 0.000012095430 1.4023330 0.1608158453402
## [2,] 0.0131376894 0.008196721 0.000017275154 1.1887789 0.2345266842014
## [3,] 0.0198717287 0.008196721 0.000017280718 2.8085119 0.0049771040414
## [4,] 0.0254654129 0.008196721 0.000014195850 4.5833048 0.0000045768422
## [5,] 0.0131335947 0.008196721 0.000014252849 1.3076790 0.1909822112521
## [6,] 0.0116851451 0.008196721 0.000021696631 0.7489160 0.4539078643927
## [7,] 0.0122618151 0.008196721 0.000014102282 1.0824950 0.2790326542547
## [8,] 0.0224384441 0.008196721 0.000013789193 3.8352450 0.0001254392000
## [9,] 0.0117686221 0.008196721 0.000017193705 0.8614194 0.3890070965404
## [10,] 0.0112492970 0.008196721 0.000017177361 0.7365262 0.4614105042275
## [11,] 0.0132424884 0.008196721 0.000044001104 0.7606685 0.4468551106494
## [12,] 0.0145077804 0.008196721 0.000012052795 1.8178513 0.0690868588962
## [13,] 0.0129126950 0.008196721 0.000016803228 1.1504692 0.2499506550092
## [14,] 0.0105583458 0.008196721 0.000029144813 0.4374518 0.6617837272740
## [15,] 0.0075028538 0.008196721 0.000021692283 -0.1489786 0.8815705047423
## [16,] 0.0137832521 0.008196721 0.000021714528 1.1988559 0.2305839685432
## [17,] 0.0258991127 0.008196721 0.000016451367 4.3644643 0.0000127434661
## [18,] 0.0164563661 0.008196721 0.000010401823 2.5609826 0.0104376569128
## [19,] 0.0184020563 0.008196721 0.000021570527 2.1973381 0.0279963073197
## [20,] 0.0112619816 0.008196721 0.000021098419 0.6673327 0.5045596417553
## [21,] 0.0155208874 0.008196721 0.000010411325 2.2698919 0.0232141412705
## [22,] 0.0157138147 0.008196721 0.000021474329 1.6221467 0.1047719322864
## [23,] 0.0157310831 0.008196721 0.000009142877 2.4917532 0.0127114316274
## [24,] 0.0164250753 0.008196721 0.000014280890 2.1773855 0.0294518151042
## [25,] 0.0236013952 0.008196721 0.000017215712 3.7127013 0.0002050588576
## [26,] 0.0121121752 0.008196721 0.000014281151 1.0360972 0.3001567987920
## [27,] 0.0277156464 0.008196721 0.000021743525 4.1859200 0.0000284013372
## [28,] 0.0165185166 0.008196721 0.000021596251 1.7907207 0.0733381308427
## [29,] 0.0143140912 0.008196721 0.000014256185 1.6201793 0.1051937713384
## [30,] 0.0296224528 0.008196721 0.000017207605 5.1650608 0.0000002403607
## [31,] 0.0116833807 0.008196721 0.000008273152 1.2121993 0.2254360931278
## [32,] 0.0012327248 0.008196721 0.000016876258 -1.6951981 0.0900378429644
## [33,] 0.0023662738 0.008196721 0.000021110840 -1.2689635 0.2044540765636
## [34,] 0.0031890616 0.008196721 0.000028338971 -0.9406819 0.3468679126656
## [35,] 0.0020310182 0.008196721 0.000028502299 -1.1548954 0.2481332545047
## [36,] 0.0026475162 0.008196721 0.000013924639 -1.4870953 0.1369896382763
## [37,] 0.0045063055 0.008196721 0.000010223325 -1.1541951 0.2484202099147
## [38,] 0.0023293255 0.008196721 0.000011768265 -1.7103664 0.0871981361078
## [39,] 0.0032337174 0.008196721 0.000014112443 -1.3211238 0.1864600830176
## [40,] 0.0051867937 0.008196721 0.000043102444 -0.4584639 0.6466192282599
## [41,] 0.0041639600 0.008196721 0.000011969961 -1.1656178 0.2437690159198
## [42,] 0.0197051854 0.008196721 0.000016815704 2.8064663 0.0050088148951
## [43,] 0.0018159417 0.008196721 0.000080644318 -0.7105373 0.4773710307771
## [44,] 0.0067350428 0.008196721 0.000028338971 -0.2745743 0.7836433366700
## [45,] 0.0048985030 0.008196721 0.000021184080 -0.7165964 0.4736231956126
## [46,] 0.0064156543 0.008196721 0.000016843871 -0.4339696 0.6643105164214
## [47,] 0.0043100453 0.008196721 0.000028838598 -0.7237545 0.4692164739028
## [48,] 0.0184008808 0.008196721 0.000017111076 2.4668263 0.0136316451239
## [49,] 0.0069255486 0.008196721 0.000028338971 -0.2387880 0.8112699620346
## [50,] 0.0044812274 0.008196721 0.000014214864 -0.9854741 0.3243912701523
## [51,] 0.0038226744 0.008196721 0.000043203523 -0.6654629 0.5057544601694
## [52,] 0.0074482126 0.008196721 0.000021266752 -0.1623104 0.8710614460200
## [53,] 0.0042997009 0.008196721 0.000021423131 -0.8419595 0.3998106071470
## [54,] 0.0138956281 0.008196721 0.000008064399 2.0068068 0.0447702424843
## [55,] 0.0044034363 0.008196721 0.000017220448 -0.9140990 0.3606648056928
## [56,] 0.0031810998 0.008196721 0.000021588239 -1.0794838 0.2803720900510
## [57,] 0.0067696966 0.008196721 0.000021240273 -0.3096360 0.7568377846068
## [58,] 0.0010634376 0.008196721 0.000021158584 -1.5507657 0.1209578337871
## [59,] 0.0052491254 0.008196721 0.000011797098 -0.8581839 0.3907909232812
## [60,] 0.0006761049 0.008196721 0.000017077047 -1.8198980 0.0687745340556
## [61,] 0.0050436403 0.008196721 0.000021075621 -0.6868232 0.4921941593385
## [62,] 0.0039620445 0.008196721 0.000017004133 -1.0269352 0.3044509877977
## [63,] 0.0020426300 0.008196721 0.000021411972 -1.3299509 0.1835344562680
## [64,] 0.0025742454 0.008196721 0.000013878565 -1.5092296 0.1312401097096
## [65,] 0.0017896203 0.008196721 0.000021653774 -1.3768758 0.1685506406435
## [66,] 0.0038725467 0.008196721 0.000021221448 -0.9386760 0.3478971116224
## [67,] 0.0031106763 0.008196721 0.000021286471 -1.1023725 0.2702997645785
## [68,] 0.0047404194 0.008196721 0.000010370453 -1.0732795 0.2831457464410
## [69,] 0.0039825165 0.008196721 0.000011834235 -1.2250267 0.2205651680142
## [70,] 0.0042307753 0.008196721 0.000011766868 -1.1561558 0.2476174385511
## [71,] 0.0042594861 0.008196721 0.000021140432 -0.8563168 0.3918225963041
## [72,] 0.0029141293 0.008196721 0.000029235427 -0.9769950 0.3285716145660
## [73,] 0.0038183277 0.008196721 0.000013833154 -1.1772107 0.2391114269497
## [74,] 0.0030737402 0.008196721 0.000021680940 -1.1002306 0.2712316804429
## [75,] 0.0035858303 0.008196721 0.000021718297 -0.9894001 0.3224674024121
## [76,] 0.0020883682 0.008196721 0.000008096533 -2.1467160 0.0318159016602
## [77,] 0.0053438992 0.008196721 0.000011786792 -0.8309538 0.4059997323292
## [78,] 0.0032662067 0.008196721 0.000021219065 -1.0703582 0.2844581187779
## [79,] 0.0037040044 0.008196721 0.000016827197 -1.0952245 0.2734182764771
## [80,] 0.0042481448 0.008196721 0.000011943137 -1.1425662 0.2532187625257
## [81,] 0.0072948707 0.008196721 0.000013812271 -0.2426622 0.8082671268726
## [82,] 0.0040181164 0.008196721 0.000017185811 -1.0079669 0.3134703310049
## [83,] 0.0014998196 0.008196721 0.000028415733 -1.2563033 0.2090060386896
## [84,] 0.0025203742 0.008196721 0.000021494883 -1.2243382 0.2208246627724
## [85,] 0.0026611627 0.008196721 0.000021321038 -1.1988291 0.2305943995008
## [86,] 0.0015638975 0.008196721 0.000017030100 -1.6072738 0.1079943173271
## [87,] 0.0016144767 0.008196721 0.000013834966 -1.7696402 0.0767870976051
## [88,] 0.0016136113 0.008196721 0.000013843756 -1.7693109 0.0768420109934
## [89,] 0.0026055040 0.008196721 0.000021075621 -1.2179128 0.2232571261932
## [90,] 0.0014914094 0.008196721 0.000028338971 -1.2595835 0.2078196617412
## [91,] 0.0017337789 0.008196721 0.000011876794 -1.8753429 0.0607455715388
## [92,] 0.0020117888 0.008196721 0.000011737029 -1.8053270 0.0710235207650
## [93,] 0.0042719028 0.008196721 0.000011737029 -1.1456198 0.2519525106859
## [94,] 0.0022778257 0.008196721 0.000029088048 -1.0974466 0.2724462122745
## [95,] 0.0028897615 0.008196721 0.000007250911 -1.9708316 0.0487431412103
## [96,] 0.0036602316 0.008196721 0.000011737029 -1.3241612 0.1854495224579
## [97,] 0.0029025136 0.008196721 0.000043081791 -0.8065924 0.4199013806112
## [98,] 0.0016762825 0.008196721 0.000021134273 -1.4183492 0.1560888325270
## [99,] 0.0003633390 0.008196721 0.000021188811 -1.7017511 0.0888020455516
## [100,] 0.0015271897 0.008196721 0.000011737029 -1.9467772 0.0515614513102
## [101,] 0.0033555751 0.008196721 0.000011782551 -1.4103553 0.1584348038729
## [102,] 0.0049479007 0.008196721 0.000021554890 -0.6997660 0.4840734551855
## [103,] 0.0064601012 0.008196721 0.000016819733 -0.4234433 0.6719718542450
## [104,] 0.0028282658 0.008196721 0.000028351628 -1.0082317 0.3133432364167
## [105,] 0.0103611343 0.008196721 0.000021075621 0.4714655 0.6373083277561
## [106,] 0.0040799227 0.008196721 0.000042940382 -0.6282416 0.5298457064391
## [107,] 0.0109014931 0.008196721 0.000013945530 0.7242910 0.4688870768594
## [108,] 0.0076859990 0.008196721 0.000010324279 -0.1589480 0.8737098788773
## [109,] 0.0154332467 0.008196721 0.000014022914 1.9324620 0.0533025136999
## [110,] 0.0095433779 0.008196721 0.000011876051 0.3907696 0.6959675118241
## [111,] 0.0140673980 0.008196721 0.000021460724 1.2672607 0.2050620754513
## [112,] 0.0106258057 0.008196721 0.000017277914 0.5843822 0.5589632305307
## [113,] 0.0138807140 0.008196721 0.000028968221 1.0560697 0.2909363503425
## [114,] 0.0088708079 0.008196721 0.000009265819 0.2214490 0.8247428182119
## [115,] 0.0070690985 0.008196721 0.000011936886 -0.3263761 0.7441398363457
## [116,] 0.0048261719 0.008196721 0.000014236485 -0.8933041 0.3716943550869
## [117,] 0.0109175072 0.008196721 0.000017161286 0.6567793 0.5113228454889
## [118,] 0.0063865137 0.008196721 0.000017213113 -0.4363135 0.6626092724731
## [119,] 0.0061135242 0.008196721 0.000087727547 -0.2224141 0.8239915098185
## [120,] 0.0099108733 0.008196721 0.000021075621 0.3733870 0.7088604661011
## [121,] 0.0165514569 0.008196721 0.000017205836 2.0141640 0.0439923206612
## [122,] 0.0133945617 0.008196721 0.000021600340 1.1183884 0.2634011543756
## [123,] 0.0132932269 0.008196721 0.000017282052 1.2259559 0.2202152756321
## attr(,"cluster")
## [1] High High High High High High Low High High High High High Low High High
## [16] High High Low High Low Low Low High High High High High High High High
## [31] High Low Low Low Low Low Low Low Low Low Low Low High Low Low
## [46] Low High High Low High Low Low Low Low High High Low Low Low Low
## [61] Low Low Low Low High Low Low Low Low Low Low High Low High High
## [76] Low Low Low Low Low Low High Low Low Low Low Low Low Low Low
## [91] Low Low Low High Low Low Low Low Low Low Low Low Low Low Low
## [106] Low Low Low Low Low Low High Low High High High High High Low Low
## [121] High High High
## Levels: Low High
## attr(,"gstari")
## [1] FALSE
## attr(,"call")
## localG(x = ssdata$NOR_TXCURR, listw = sscastlistw, zero.policy = TRUE)
## attr(,"class")
## [1] "localG"Gimap <- cbind(ssdata, as.matrix(Gis), as.matrix(attributes(Gis)$internals), as.matrix(attributes(Gis)$cluster))
Gimap <- Gimap %>% rename('prz'= 'Pr.z....E.Gi..')
Gimap <- Gimap %>% rename_at('as.matrix.attributes.Gis..cluster.', ~'cluster')
tm_shape(Gimap) + tm_fill(col = "as.matrix.Gis.", palette = "-viridis", midpoint = NA, title = "G*") + tm_legend(position = c("centre", "top"))
tm_shape(Gimap) + tm_fill(col = "Gi", palette = "-viridis", title = "Local G Observed Statistic")
tm_shape(Gimap) + tm_fill(col = "prz", style ="quantile", palette = "-viridis", title = "Local G P-Value") + tm_layout(legend.position = c(0.45, 0.65))
Aspatial multivariate analysis using Poisson distribution
Poisson regression
poisson_output <-glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK,
data = rschdata, family = poisson)## Warning in dpois(y, mu, log = TRUE): non-integer x = 8481.401340## Warning in dpois(y, mu, log = TRUE): non-integer x = 5871.310184## Warning in dpois(y, mu, log = TRUE): non-integer x = 4959.670205## Warning in dpois(y, mu, log = TRUE): non-integer x = 5471.895840## Warning in dpois(y, mu, log = TRUE): non-integer x = 8420.317920## Warning in dpois(y, mu, log = TRUE): non-integer x = 2038.665138## Warning in dpois(y, mu, log = TRUE): non-integer x = 7014.040846## Warning in dpois(y, mu, log = TRUE): non-integer x = 10895.600778## Warning in dpois(y, mu, log = TRUE): non-integer x = 2953.330215## Warning in dpois(y, mu, log = TRUE): non-integer x = 2460.296866## Warning in dpois(y, mu, log = TRUE): non-integer x = 4479.285188## Warning in dpois(y, mu, log = TRUE): non-integer x = 3502.882481## Warning in dpois(y, mu, log = TRUE): non-integer x = 2691.169494## Warning in dpois(y, mu, log = TRUE): non-integer x = 4673.598586## Warning in dpois(y, mu, log = TRUE): non-integer x = 10005.176030## Warning in dpois(y, mu, log = TRUE): non-integer x = 8248.138284## Warning in dpois(y, mu, log = TRUE): non-integer x = 3740.698556## Warning in dpois(y, mu, log = TRUE): non-integer x = 5035.542691## Warning in dpois(y, mu, log = TRUE): non-integer x = 3817.352694## Warning in dpois(y, mu, log = TRUE): non-integer x = 3264.867055## Warning in dpois(y, mu, log = TRUE): non-integer x = 6173.078593## Warning in dpois(y, mu, log = TRUE): non-integer x = 6754.753217## Warning in dpois(y, mu, log = TRUE): non-integer x = 3358.161900## Warning in dpois(y, mu, log = TRUE): non-integer x = 2865.751974## Warning in dpois(y, mu, log = TRUE): non-integer x = 3451.472782## Warning in dpois(y, mu, log = TRUE): non-integer x = 1818.172777## Warning in dpois(y, mu, log = TRUE): non-integer x = 2166.468698## Warning in dpois(y, mu, log = TRUE): non-integer x = 1204.255817## Warning in dpois(y, mu, log = TRUE): non-integer x = 7128.989005## Warning in dpois(y, mu, log = TRUE): non-integer x = 3038.513644## Warning in dpois(y, mu, log = TRUE): non-integer x = 3410.668567## Warning in dpois(y, mu, log = TRUE): non-integer x = 13298.755187## Warning in dpois(y, mu, log = TRUE): non-integer x = 2854.128317## Warning in dpois(y, mu, log = TRUE): non-integer x = 3165.654120## Warning in dpois(y, mu, log = TRUE): non-integer x = 2378.993253## Warning in dpois(y, mu, log = TRUE): non-integer x = 4920.580700## Warning in dpois(y, mu, log = TRUE): non-integer x = 3083.709932## Warning in dpois(y, mu, log = TRUE): non-integer x = 4380.800637## Warning in dpois(y, mu, log = TRUE): non-integer x = 1948.638019## Warning in dpois(y, mu, log = TRUE): non-integer x = 3461.228178## Warning in dpois(y, mu, log = TRUE): non-integer x = 2038.418543## Warning in dpois(y, mu, log = TRUE): non-integer x = 989.605252## Warning in dpois(y, mu, log = TRUE): non-integer x = 2097.851365## Warning in dpois(y, mu, log = TRUE): non-integer x = 1317.327670## Warning in dpois(y, mu, log = TRUE): non-integer x = 1635.274435## Warning in dpois(y, mu, log = TRUE): non-integer x = 1512.785015## Warning in dpois(y, mu, log = TRUE): non-integer x = 3176.829154## Warning in dpois(y, mu, log = TRUE): non-integer x = 6305.475941## Warning in dpois(y, mu, log = TRUE): non-integer x = 3197.854925## Warning in dpois(y, mu, log = TRUE): non-integer x = 2766.486275## Warning in dpois(y, mu, log = TRUE): non-integer x = 3229.232904## Warning in dpois(y, mu, log = TRUE): non-integer x = 1389.514188## Warning in dpois(y, mu, log = TRUE): non-integer x = 788.023557## Warning in dpois(y, mu, log = TRUE): non-integer x = 2271.334457## Warning in dpois(y, mu, log = TRUE): non-integer x = 2156.529404## Warning in dpois(y, mu, log = TRUE): non-integer x = 993.581464## Warning in dpois(y, mu, log = TRUE): non-integer x = 1706.962657## Warning in dpois(y, mu, log = TRUE): non-integer x = 2432.547536## Warning in dpois(y, mu, log = TRUE): non-integer x = 1444.338618## Warning in dpois(y, mu, log = TRUE): non-integer x = 3187.785293## Warning in dpois(y, mu, log = TRUE): non-integer x = 2264.150943## Warning in dpois(y, mu, log = TRUE): non-integer x = 1812.692485## Warning in dpois(y, mu, log = TRUE): non-integer x = 463.902066## Warning in dpois(y, mu, log = TRUE): non-integer x = 1557.858493## Warning in dpois(y, mu, log = TRUE): non-integer x = 2708.026111## Warning in dpois(y, mu, log = TRUE): non-integer x = 819.818382## Warning in dpois(y, mu, log = TRUE): non-integer x = 2699.916909## Warning in dpois(y, mu, log = TRUE): non-integer x = 1681.717813## Warning in dpois(y, mu, log = TRUE): non-integer x = 483.908382## Warning in dpois(y, mu, log = TRUE): non-integer x = 987.692308## Warning in dpois(y, mu, log = TRUE): non-integer x = 2010.017072## Warning in dpois(y, mu, log = TRUE): non-integer x = 1676.891283## Warning in dpois(y, mu, log = TRUE): non-integer x = 630.986256## Warning in dpois(y, mu, log = TRUE): non-integer x = 1274.897249## Warning in dpois(y, mu, log = TRUE): non-integer x = 627.822492## Warning in dpois(y, mu, log = TRUE): non-integer x = 301.822941## Warning in dpois(y, mu, log = TRUE): non-integer x = 521.197867## Warning in dpois(y, mu, log = TRUE): non-integer x = 736.868232## Warning in dpois(y, mu, log = TRUE): non-integer x = 969.112947## Warning in dpois(y, mu, log = TRUE): non-integer x = 774.004910## Warning in dpois(y, mu, log = TRUE): non-integer x = 633.764206## Warning in dpois(y, mu, log = TRUE): non-integer x = 668.274128## Warning in dpois(y, mu, log = TRUE): non-integer x = 605.577233## Warning in dpois(y, mu, log = TRUE): non-integer x = 655.602657## Warning in dpois(y, mu, log = TRUE): non-integer x = 402.460159## Warning in dpois(y, mu, log = TRUE): non-integer x = 435.482057## Warning in dpois(y, mu, log = TRUE): non-integer x = 208.910468## Warning in dpois(y, mu, log = TRUE): non-integer x = 404.490875## Warning in dpois(y, mu, log = TRUE): non-integer x = 330.872183## Warning in dpois(y, mu, log = TRUE): non-integer x = 456.819556## Warning in dpois(y, mu, log = TRUE): non-integer x = 551.828945## Warning in dpois(y, mu, log = TRUE): non-integer x = 376.207687## Warning in dpois(y, mu, log = TRUE): non-integer x = 327.422414## Warning in dpois(y, mu, log = TRUE): non-integer x = 542.296508## Warning in dpois(y, mu, log = TRUE): non-integer x = 395.957829## Warning in dpois(y, mu, log = TRUE): non-integer x = 205.185600## Warning in dpois(y, mu, log = TRUE): non-integer x = 402.178212## Warning in dpois(y, mu, log = TRUE): non-integer x = 127.888291## Warning in dpois(y, mu, log = TRUE): non-integer x = 455.271689## Warning in dpois(y, mu, log = TRUE): non-integer x = 237.828891## Warning in dpois(y, mu, log = TRUE): non-integer x = 422.899080## Warning in dpois(y, mu, log = TRUE): non-integer x = 214.751508## Warning in dpois(y, mu, log = TRUE): non-integer x = 306.464857## Warning in dpois(y, mu, log = TRUE): non-integer x = 195.831658## Warning in dpois(y, mu, log = TRUE): non-integer x = 125.718017## Warning in dpois(y, mu, log = TRUE): non-integer x = 133.459532## Warning in dpois(y, mu, log = TRUE): non-integer x = 175.237327## Warning in dpois(y, mu, log = TRUE): non-integer x = 82.428390## Warning in dpois(y, mu, log = TRUE): non-integer x = 115.570751## Warning in dpois(y, mu, log = TRUE): non-integer x = 33.879776print(summary(poisson_output))##
## Call:
## glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV +
## HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, family = poisson,
## data = rschdata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -97.790 -17.152 -0.886 17.732 83.428
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.22781639879 0.00628298986 1309.54 <0.0000000000000002 ***
## POP_DEN 0.00006018465 0.00000219290 27.45 <0.0000000000000002 ***
## PHF_COV 0.01568730882 0.00016649608 94.22 <0.0000000000000002 ***
## AHF_COV -0.00618443118 0.00003826263 -161.63 <0.0000000000000002 ***
## POP_POV -0.00000317098 0.00000005951 -53.28 <0.0000000000000002 ***
## HEALTH_DEP 0.00000701215 0.00000018534 37.83 <0.0000000000000002 ***
## EDUC_DEP -0.00003776223 0.00000031499 -119.89 <0.0000000000000002 ***
## WORK_DEP 0.00002833237 0.00000026063 108.71 <0.0000000000000002 ***
## SEC_SHOCK -0.00006747833 0.00000042840 -157.51 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 304081 on 122 degrees of freedom
## Residual deviance: 136424 on 114 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 5vif_values <- vif(poisson_output)
vif_values## POP_DEN PHF_COV AHF_COV POP_POV HEALTH_DEP EDUC_DEP WORK_DEP
## 1.614147 1.360960 1.072502 8.350692 22.275535 3.684154 14.762289
## SEC_SHOCK
## 3.907127barplot(vif_values, main = "VIF Values", horiz = TRUE, col = "black", cex.names = 0.5)+ abline(v = 5, lwd = 3, lty = 2)
## numeric(0)plot(poisson_output)
poisson_output1 <-glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + EDUC_DEP + WORK_DEP + SEC_SHOCK,
data = rschdata, family = poisson)## Warning in dpois(y, mu, log = TRUE): non-integer x = 8481.401340## Warning in dpois(y, mu, log = TRUE): non-integer x = 5871.310184## Warning in dpois(y, mu, log = TRUE): non-integer x = 4959.670205## Warning in dpois(y, mu, log = TRUE): non-integer x = 5471.895840## Warning in dpois(y, mu, log = TRUE): non-integer x = 8420.317920## Warning in dpois(y, mu, log = TRUE): non-integer x = 2038.665138## Warning in dpois(y, mu, log = TRUE): non-integer x = 7014.040846## Warning in dpois(y, mu, log = TRUE): non-integer x = 10895.600778## Warning in dpois(y, mu, log = TRUE): non-integer x = 2953.330215## Warning in dpois(y, mu, log = TRUE): non-integer x = 2460.296866## Warning in dpois(y, mu, log = TRUE): non-integer x = 4479.285188## Warning in dpois(y, mu, log = TRUE): non-integer x = 3502.882481## Warning in dpois(y, mu, log = TRUE): non-integer x = 2691.169494## Warning in dpois(y, mu, log = TRUE): non-integer x = 4673.598586## Warning in dpois(y, mu, log = TRUE): non-integer x = 10005.176030## Warning in dpois(y, mu, log = TRUE): non-integer x = 8248.138284## Warning in dpois(y, mu, log = TRUE): non-integer x = 3740.698556## Warning in dpois(y, mu, log = TRUE): non-integer x = 5035.542691## Warning in dpois(y, mu, log = TRUE): non-integer x = 3817.352694## Warning in dpois(y, mu, log = TRUE): non-integer x = 3264.867055## Warning in dpois(y, mu, log = TRUE): non-integer x = 6173.078593## Warning in dpois(y, mu, log = TRUE): non-integer x = 6754.753217## Warning in dpois(y, mu, log = TRUE): non-integer x = 3358.161900## Warning in dpois(y, mu, log = TRUE): non-integer x = 2865.751974## Warning in dpois(y, mu, log = TRUE): non-integer x = 3451.472782## Warning in dpois(y, mu, log = TRUE): non-integer x = 1818.172777## Warning in dpois(y, mu, log = TRUE): non-integer x = 2166.468698## Warning in dpois(y, mu, log = TRUE): non-integer x = 1204.255817## Warning in dpois(y, mu, log = TRUE): non-integer x = 7128.989005## Warning in dpois(y, mu, log = TRUE): non-integer x = 3038.513644## Warning in dpois(y, mu, log = TRUE): non-integer x = 3410.668567## Warning in dpois(y, mu, log = TRUE): non-integer x = 13298.755187## Warning in dpois(y, mu, log = TRUE): non-integer x = 2854.128317## Warning in dpois(y, mu, log = TRUE): non-integer x = 3165.654120## Warning in dpois(y, mu, log = TRUE): non-integer x = 2378.993253## Warning in dpois(y, mu, log = TRUE): non-integer x = 4920.580700## Warning in dpois(y, mu, log = TRUE): non-integer x = 3083.709932## Warning in dpois(y, mu, log = TRUE): non-integer x = 4380.800637## Warning in dpois(y, mu, log = TRUE): non-integer x = 1948.638019## Warning in dpois(y, mu, log = TRUE): non-integer x = 3461.228178## Warning in dpois(y, mu, log = TRUE): non-integer x = 2038.418543## Warning in dpois(y, mu, log = TRUE): non-integer x = 989.605252## Warning in dpois(y, mu, log = TRUE): non-integer x = 2097.851365## Warning in dpois(y, mu, log = TRUE): non-integer x = 1317.327670## Warning in dpois(y, mu, log = TRUE): non-integer x = 1635.274435## Warning in dpois(y, mu, log = TRUE): non-integer x = 1512.785015## Warning in dpois(y, mu, log = TRUE): non-integer x = 3176.829154## Warning in dpois(y, mu, log = TRUE): non-integer x = 6305.475941## Warning in dpois(y, mu, log = TRUE): non-integer x = 3197.854925## Warning in dpois(y, mu, log = TRUE): non-integer x = 2766.486275## Warning in dpois(y, mu, log = TRUE): non-integer x = 3229.232904## Warning in dpois(y, mu, log = TRUE): non-integer x = 1389.514188## Warning in dpois(y, mu, log = TRUE): non-integer x = 788.023557## Warning in dpois(y, mu, log = TRUE): non-integer x = 2271.334457## Warning in dpois(y, mu, log = TRUE): non-integer x = 2156.529404## Warning in dpois(y, mu, log = TRUE): non-integer x = 993.581464## Warning in dpois(y, mu, log = TRUE): non-integer x = 1706.962657## Warning in dpois(y, mu, log = TRUE): non-integer x = 2432.547536## Warning in dpois(y, mu, log = TRUE): non-integer x = 1444.338618## Warning in dpois(y, mu, log = TRUE): non-integer x = 3187.785293## Warning in dpois(y, mu, log = TRUE): non-integer x = 2264.150943## Warning in dpois(y, mu, log = TRUE): non-integer x = 1812.692485## Warning in dpois(y, mu, log = TRUE): non-integer x = 463.902066## Warning in dpois(y, mu, log = TRUE): non-integer x = 1557.858493## Warning in dpois(y, mu, log = TRUE): non-integer x = 2708.026111## Warning in dpois(y, mu, log = TRUE): non-integer x = 819.818382## Warning in dpois(y, mu, log = TRUE): non-integer x = 2699.916909## Warning in dpois(y, mu, log = TRUE): non-integer x = 1681.717813## Warning in dpois(y, mu, log = TRUE): non-integer x = 483.908382## Warning in dpois(y, mu, log = TRUE): non-integer x = 987.692308## Warning in dpois(y, mu, log = TRUE): non-integer x = 2010.017072## Warning in dpois(y, mu, log = TRUE): non-integer x = 1676.891283## Warning in dpois(y, mu, log = TRUE): non-integer x = 630.986256## Warning in dpois(y, mu, log = TRUE): non-integer x = 1274.897249## Warning in dpois(y, mu, log = TRUE): non-integer x = 627.822492## Warning in dpois(y, mu, log = TRUE): non-integer x = 301.822941## Warning in dpois(y, mu, log = TRUE): non-integer x = 521.197867## Warning in dpois(y, mu, log = TRUE): non-integer x = 736.868232## Warning in dpois(y, mu, log = TRUE): non-integer x = 969.112947## Warning in dpois(y, mu, log = TRUE): non-integer x = 774.004910## Warning in dpois(y, mu, log = TRUE): non-integer x = 633.764206## Warning in dpois(y, mu, log = TRUE): non-integer x = 668.274128## Warning in dpois(y, mu, log = TRUE): non-integer x = 605.577233## Warning in dpois(y, mu, log = TRUE): non-integer x = 655.602657## Warning in dpois(y, mu, log = TRUE): non-integer x = 402.460159## Warning in dpois(y, mu, log = TRUE): non-integer x = 435.482057## Warning in dpois(y, mu, log = TRUE): non-integer x = 208.910468## Warning in dpois(y, mu, log = TRUE): non-integer x = 404.490875## Warning in dpois(y, mu, log = TRUE): non-integer x = 330.872183## Warning in dpois(y, mu, log = TRUE): non-integer x = 456.819556## Warning in dpois(y, mu, log = TRUE): non-integer x = 551.828945## Warning in dpois(y, mu, log = TRUE): non-integer x = 376.207687## Warning in dpois(y, mu, log = TRUE): non-integer x = 327.422414## Warning in dpois(y, mu, log = TRUE): non-integer x = 542.296508## Warning in dpois(y, mu, log = TRUE): non-integer x = 395.957829## Warning in dpois(y, mu, log = TRUE): non-integer x = 205.185600## Warning in dpois(y, mu, log = TRUE): non-integer x = 402.178212## Warning in dpois(y, mu, log = TRUE): non-integer x = 127.888291## Warning in dpois(y, mu, log = TRUE): non-integer x = 455.271689## Warning in dpois(y, mu, log = TRUE): non-integer x = 237.828891## Warning in dpois(y, mu, log = TRUE): non-integer x = 422.899080## Warning in dpois(y, mu, log = TRUE): non-integer x = 214.751508## Warning in dpois(y, mu, log = TRUE): non-integer x = 306.464857## Warning in dpois(y, mu, log = TRUE): non-integer x = 195.831658## Warning in dpois(y, mu, log = TRUE): non-integer x = 125.718017## Warning in dpois(y, mu, log = TRUE): non-integer x = 133.459532## Warning in dpois(y, mu, log = TRUE): non-integer x = 175.237327## Warning in dpois(y, mu, log = TRUE): non-integer x = 82.428390## Warning in dpois(y, mu, log = TRUE): non-integer x = 115.570751## Warning in dpois(y, mu, log = TRUE): non-integer x = 33.879776print(summary(poisson_output1))##
## Call:
## glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV +
## EDUC_DEP + WORK_DEP + SEC_SHOCK, family = poisson, data = rschdata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -98.621 -19.121 0.351 18.421 81.594
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.31858168385 0.00583431141 1425.80 <0.0000000000000002 ***
## POP_DEN 0.00008990058 0.00000203674 44.14 <0.0000000000000002 ***
## PHF_COV 0.01400305796 0.00016128267 86.82 <0.0000000000000002 ***
## AHF_COV -0.00626078615 0.00003784892 -165.41 <0.0000000000000002 ***
## POP_POV -0.00000216883 0.00000005202 -41.69 <0.0000000000000002 ***
## EDUC_DEP -0.00003102203 0.00000025882 -119.86 <0.0000000000000002 ***
## WORK_DEP 0.00003336394 0.00000023380 142.70 <0.0000000000000002 ***
## SEC_SHOCK -0.00006514925 0.00000042239 -154.24 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 304081 on 122 degrees of freedom
## Residual deviance: 137907 on 115 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 5vif_values1 <- vif(poisson_output1)
vif_values1## POP_DEN PHF_COV AHF_COV POP_POV EDUC_DEP WORK_DEP SEC_SHOCK
## 1.376051 1.254362 1.062385 6.972575 2.511460 12.106627 3.837199library(MASS)## Warning: package 'MASS' was built under R version 4.2.3##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectgnb <- glm(NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, family = negative.binomial(2), data = rschdata)
print(summary(gnb))##
## Call:
## glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV +
## HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, family = negative.binomial(2),
## data = rschdata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.5655 -1.0226 -0.1408 0.5074 2.2560
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.055704222 0.248862370 32.370 < 0.0000000000000002 ***
## POP_DEN 0.000213907 0.000099261 2.155 0.03327 *
## PHF_COV 0.006443845 0.007806149 0.825 0.41082
## AHF_COV -0.003067233 0.000383892 -7.990 0.0000000000012 ***
## POP_POV -0.000002114 0.000001927 -1.097 0.27508
## HEALTH_DEP 0.000002879 0.000006242 0.461 0.64544
## EDUC_DEP -0.000027998 0.000010417 -2.688 0.00827 **
## WORK_DEP 0.000018453 0.000006549 2.818 0.00570 **
## SEC_SHOCK -0.000040491 0.000012146 -3.334 0.00116 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2) family taken to be 1.286706)
##
## Null deviance: 628.15 on 122 degrees of freedom
## Residual deviance: 466.13 on 114 degrees of freedom
## AIC: 2204.7
##
## Number of Fisher Scoring iterations: 24pchisq(466.25, 114, lower.tail = FALSE)## [1] 0.00000000000000000000000000000000000000000004043444lm.morantest(gnb, sscastlistw, alternative = "two.sided")##
## Global Moran I for regression residuals
##
## data:
## model: glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV
## + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, family =
## negative.binomial(2), data = rschdata)
## weights: sscastlistw
##
## Moran I statistic standard deviate = 9.6301, p-value <
## 0.00000000000000022
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I Expectation Variance
## 0.559071077 -0.008799756 0.003477290lm.LMtests(gnb, sscastlistw, test = c("LMerr", "LMlag", "RLMerr", "RLMlag", "SARMA"))##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV
## + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, family =
## negative.binomial(2), data = rschdata)
## weights: sscastlistw
##
## LMerr = 90.124, df = 1, p-value < 0.00000000000000022
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV
## + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, family =
## negative.binomial(2), data = rschdata)
## weights: sscastlistw
##
## LMlag = 2.9942, df = 1, p-value = 0.08356
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV
## + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, family =
## negative.binomial(2), data = rschdata)
## weights: sscastlistw
##
## RLMerr = 90.099, df = 1, p-value < 0.00000000000000022
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV
## + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, family =
## negative.binomial(2), data = rschdata)
## weights: sscastlistw
##
## RLMlag = 2.9692, df = 1, p-value = 0.08487
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV
## + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, family =
## negative.binomial(2), data = rschdata)
## weights: sscastlistw
##
## SARMA = 93.093, df = 2, p-value < 0.00000000000000022Geographical weighted poisson regression
ssdataspdf <- sf:::as_Spatial(ssdata)
DM<-gw.dist(dp.locat=coordinates(ssdataspdf))
ssdataspdf$NOR_TXCURR <- round(ssdataspdf$NOR_TXCURR, digits = 0)
bw <- bw.ggwr(NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, data = ssdataspdf, family = "poisson", approach = "AICc", kernel = "exponential", adaptive = TRUE, dMat = DM )## Iteration Log-Likelihood(With bandwidth: 83 )
## =========================
## 0 -2.325e+05
## 1 -2.779e+09
## 2 -1.023e+09
## 3 -3.762e+08
## 4 -1.385e+08
## 5 -5.1e+07
## 6 -1.882e+07
## 7 -6.974e+06
## 8 -2.603e+06
## 9 -9.892e+05
## 10 -3.912e+05
## 11 -1.701e+05
## 12 -9.073e+04
## 13 -6.528e+04
## 14 -5.9e+04
## 15 -5.761e+04
## 16 -5.728e+04
## 17 -5.727e+04
## 18 -5.727e+04
## 19 -5.727e+04
## Adaptive bandwidth (number of nearest neighbours): 83 AICc value: 114566.1
## Iteration Log-Likelihood(With bandwidth: 59 )
## =========================
## 0 -1.966e+05
## 1 -2.149e+07
## 2 -7.924e+06
## 3 -2.915e+06
## 4 -1.082e+06
## 5 -4.158e+05
## 6 -1.759e+05
## 7 -9.096e+04
## 8 -6.338e+04
## 9 -5.613e+04
## 10 -5.437e+04
## 11 -5.397e+04
## 12 -5.394e+04
## 13 -5.394e+04
## 14 -5.394e+04
## Adaptive bandwidth (number of nearest neighbours): 59 AICc value: 107922.9
## Iteration Log-Likelihood(With bandwidth: 43 )
## =========================
## 0 -1.612e+05
## 1 -1.572e+06
## 2 -5.604e+05
## 3 -2.083e+05
## 4 -9.491e+04
## 5 -6.116e+04
## 6 -5.173e+04
## 7 -4.983e+04
## 8 -4.953e+04
## 9 -4.949e+04
## 10 -4.949e+04
## Adaptive bandwidth (number of nearest neighbours): 43 AICc value: 99024.94
## Iteration Log-Likelihood(With bandwidth: 34 )
## =========================
## 0 -1.412e+05
## 1 -1.484e+06
## 2 -5.596e+05
## 3 -2.169e+05
## 4 -1.004e+05
## 5 -6.217e+04
## 6 -5.001e+04
## 7 -4.689e+04
## 8 -4.631e+04
## 9 -4.619e+04
## 10 -4.618e+04
## 11 -4.618e+04
## Adaptive bandwidth (number of nearest neighbours): 34 AICc value: 92424.72
## Iteration Log-Likelihood(With bandwidth: 27 )
## =========================
## 0 -1.277e+05
## 1 -1.953e+06
## 2 -7.668e+05
## 3 -2.979e+05
## 4 -1.295e+05
## 5 -7.1e+04
## 6 -5.099e+04
## 7 -4.497e+04
## 8 -4.367e+04
## 9 -4.34e+04
## 10 -4.336e+04
## 11 -4.336e+04
## Adaptive bandwidth (number of nearest neighbours): 27 AICc value: 86784.45
## Iteration Log-Likelihood(With bandwidth: 24 )
## =========================
## 0 -1.217e+05
## 1 -2.497e+06
## 2 -1.003e+06
## 3 -3.909e+05
## 4 -1.653e+05
## 5 -8.386e+04
## 6 -5.435e+04
## 7 -4.475e+04
## 8 -4.239e+04
## 9 -4.191e+04
## 10 -4.181e+04
## 11 -4.181e+04
## 12 -4.181e+04
## 13 -4.181e+04
## Adaptive bandwidth (number of nearest neighbours): 24 AICc value: 83690.03
## Iteration Log-Likelihood(With bandwidth: 21 )
## =========================
## 0 -1.15e+05
## 1 -2.547e+06
## 2 -1.016e+06
## 3 -3.939e+05
## 4 -1.658e+05
## 5 -8.316e+04
## 6 -5.345e+04
## 7 -4.363e+04
## 8 -4.103e+04
## 9 -4.047e+04
## 10 -4.035e+04
## 11 -4.035e+04
## 12 -4.035e+04
## 13 -4.035e+04
## Adaptive bandwidth (number of nearest neighbours): 21 AICc value: 80771.03
## Iteration Log-Likelihood(With bandwidth: 20 )
## =========================
## 0 -1.14e+05
## 1 -2.556e+06
## 2 -1.018e+06
## 3 -3.937e+05
## 4 -1.649e+05
## 5 -8.232e+04
## 6 -5.282e+04
## 7 -4.313e+04
## 8 -4.057e+04
## 9 -4.001e+04
## 10 -3.989e+04
## 11 -3.988e+04
## 12 -3.988e+04
## 13 -3.988e+04
## Adaptive bandwidth (number of nearest neighbours): 20 AICc value: 79844.17
## Iteration Log-Likelihood(With bandwidth: 18 )
## =========================
## 0 -1.087e+05
## 1 -2.794e+06
## 2 -1.121e+06
## 3 -4.336e+05
## 4 -1.793e+05
## 5 -8.607e+04
## 6 -5.274e+04
## 7 -4.193e+04
## 8 -3.908e+04
## 9 -3.846e+04
## 10 -3.832e+04
## 11 -3.831e+04
## 12 -3.831e+04
## 13 -3.831e+04
## Adaptive bandwidth (number of nearest neighbours): 18 AICc value: 76710.47
## Iteration Log-Likelihood(With bandwidth: 18 )
## =========================
## 0 -1.087e+05
## 1 -2.794e+06
## 2 -1.121e+06
## 3 -4.336e+05
## 4 -1.793e+05
## 5 -8.607e+04
## 6 -5.274e+04
## 7 -4.193e+04
## 8 -3.908e+04
## 9 -3.846e+04
## 10 -3.832e+04
## 11 -3.831e+04
## 12 -3.831e+04
## 13 -3.831e+04
## Adaptive bandwidth (number of nearest neighbours): 18 AICc value: 76710.47gwr <- ggwr.basic(NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, data = ssdataspdf, bw=bw, family ="poisson", kernel = "exponential", adaptive = FALSE, dMat=DM)## Iteration Log-Likelihood
## =========================
## 0 -3.286e+05
## 1 -2.356e+11
## 2 -8.667e+10
## 3 -3.189e+10
## 4 -1.173e+10
## 5 -4.315e+09
## 6 -1.588e+09
## 7 -5.844e+08
## 8 -2.152e+08
## 9 -7.932e+07
## 10 -2.931e+07
## 11 -1.088e+07
## 12 -4.071e+06
## 13 -1.548e+06
## 14 -6.07e+05
## 15 -2.545e+05
## 16 -1.242e+05
## 17 -7.989e+04
## 18 -6.872e+04
## 19 -6.726e+04gwr## ***********************************************************************
## * Package GWmodel *
## ***********************************************************************
## Program starts at: 2025-01-11 18:20:39
## Call:
## ggwr.basic(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV +
## POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, data = ssdataspdf,
## bw = bw, family = "poisson", kernel = "exponential", adaptive = FALSE,
## dMat = DM)
##
## Dependent (y) variable: NOR_TXCURR
## Independent variables: POP_DEN PHF_COV AHF_COV POP_POV HEALTH_DEP EDUC_DEP WORK_DEP SEC_SHOCK
## Number of data points: 123
## Used family: poisson
## ***********************************************************************
## * Results of Generalized linear Regression *
## ***********************************************************************
##
## Call:
## NULL
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.000 -0.387 -0.017 0.615 279.945
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## Intercept 8.22779520411 0.00628295755 1309.54 <0.0000000000000002 ***
## POP_DEN 0.00006017078 0.00000219293 27.44 <0.0000000000000002 ***
## PHF_COV 0.01568763771 0.00016649451 94.22 <0.0000000000000002 ***
## AHF_COV -0.00618418686 0.00003826096 -161.63 <0.0000000000000002 ***
## POP_POV -0.00000317094 0.00000005951 -53.28 <0.0000000000000002 ***
## HEALTH_DEP 0.00000701185 0.00000018534 37.83 <0.0000000000000002 ***
## EDUC_DEP -0.00003776003 0.00000031498 -119.88 <0.0000000000000002 ***
## WORK_DEP 0.00002833128 0.00000026063 108.70 <0.0000000000000002 ***
## SEC_SHOCK -0.00006747452 0.00000042839 -157.51 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 304084 on 122 degrees of freedom
## Residual deviance: 136437 on 114 degrees of freedom
## AIC: 136455
##
## Number of Fisher Scoring iterations: 5
##
##
## AICc: 136456.5
## Pseudo R-square value: 0.5513183
## ***********************************************************************
## * Results of Geographically Weighted Regression *
## ***********************************************************************
##
## *********************Model calibration information*********************
## Kernel function: exponential
## Fixed bandwidth: 18
## Regression points: the same locations as observations are used.
## Distance metric: A distance matrix is specified for this model calibration.
##
## ************Summary of Generalized GWR coefficient estimates:**********
## Min. 1st Qu. Median 3rd Qu. Max.
## Intercept 8.2111747227 8.2187269022 8.2417798013 8.2560958771 8.2705
## POP_DEN 0.0000463682 0.0000481389 0.0000516193 0.0000573589 0.0001
## PHF_COV 0.0149900554 0.0152600207 0.0155711444 0.0157040793 0.0158
## AHF_COV -0.0062366651 -0.0061907384 -0.0061425612 -0.0060675648 -0.0060
## POP_POV -0.0000033274 -0.0000032653 -0.0000032028 -0.0000030079 0.0000
## HEALTH_DEP 0.0000068131 0.0000068471 0.0000069027 0.0000069570 0.0000
## EDUC_DEP -0.0000372629 -0.0000370117 -0.0000365966 -0.0000362791 0.0000
## WORK_DEP 0.0000271300 0.0000274056 0.0000282319 0.0000285182 0.0000
## SEC_SHOCK -0.0000679106 -0.0000675655 -0.0000673194 -0.0000668842 -0.0001
## ************************Diagnostic information*************************
## Number of data points: 123
## GW Deviance: 133534.9
## AIC : 133554
## AICc : 133555.7
## Pseudo R-square value: 0.5608618
##
## ***********************************************************************
## Program stops at: 2025-01-11 18:20:39bw1 <- bw.ggwr(NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + EDUC_DEP + WORK_DEP + SEC_SHOCK, data = ssdataspdf, family = "poisson", approach = "AICc", kernel = "exponential", adaptive = TRUE, dMat = DM )## Iteration Log-Likelihood(With bandwidth: 83 )
## =========================
## 0 -2.282e+05
## 1 -1.359e+10
## 2 -5e+09
## 3 -1.839e+09
## 4 -6.768e+08
## 5 -2.491e+08
## 6 -9.173e+07
## 7 -3.382e+07
## 8 -1.251e+07
## 9 -4.654e+06
## 10 -1.753e+06
## 11 -6.765e+05
## 12 -2.768e+05
## 13 -1.292e+05
## 14 -7.76e+04
## 15 -6.242e+04
## 16 -5.903e+04
## 17 -5.818e+04
## 18 -5.805e+04
## 19 -5.805e+04
## Adaptive bandwidth (number of nearest neighbours): 83 AICc value: 116124.9
## Iteration Log-Likelihood(With bandwidth: 59 )
## =========================
## 0 -1.936e+05
## 1 -1.654e+08
## 2 -6.095e+07
## 3 -2.242e+07
## 4 -8.265e+06
## 5 -3.06e+06
## 6 -1.151e+06
## 7 -4.467e+05
## 8 -1.879e+05
## 9 -9.569e+04
## 10 -6.538e+04
## 11 -5.727e+04
## 12 -5.532e+04
## 13 -5.485e+04
## 14 -5.482e+04
## 15 -5.482e+04
## 16 -5.482e+04
## Adaptive bandwidth (number of nearest neighbours): 59 AICc value: 109676.9
## Iteration Log-Likelihood(With bandwidth: 43 )
## =========================
## 0 -1.601e+05
## 1 -6.701e+06
## 2 -2.417e+06
## 3 -8.752e+05
## 4 -3.323e+05
## 5 -1.463e+05
## 6 -8.156e+04
## 7 -5.892e+04
## 8 -5.221e+04
## 9 -5.089e+04
## 10 -5.066e+04
## 11 -5.063e+04
## 12 -5.063e+04
## 13 -5.063e+04
## Adaptive bandwidth (number of nearest neighbours): 43 AICc value: 101308.6
## Iteration Log-Likelihood(With bandwidth: 34 )
## =========================
## 0 -1.4e+05
## 1 -3.842e+06
## 2 -1.412e+06
## 3 -5.254e+05
## 4 -2.153e+05
## 5 -1.067e+05
## 6 -6.705e+04
## 7 -5.288e+04
## 8 -4.858e+04
## 9 -4.767e+04
## 10 -4.75e+04
## 11 -4.748e+04
## 12 -4.748e+04
## 13 -4.748e+04
## Adaptive bandwidth (number of nearest neighbours): 34 AICc value: 95014.65
## Iteration Log-Likelihood(With bandwidth: 27 )
## =========================
## 0 -1.277e+05
## 1 -3.227e+06
## 2 -1.219e+06
## 3 -4.624e+05
## 4 -1.919e+05
## 5 -9.574e+04
## 6 -6.158e+04
## 7 -4.946e+04
## 8 -4.579e+04
## 9 -4.5e+04
## 10 -4.484e+04
## 11 -4.482e+04
## 12 -4.482e+04
## 13 -4.482e+04
## Adaptive bandwidth (number of nearest neighbours): 27 AICc value: 89705.99
## Iteration Log-Likelihood(With bandwidth: 24 )
## =========================
## 0 -1.218e+05
## 1 -2.85e+06
## 2 -1.109e+06
## 3 -4.265e+05
## 4 -1.778e+05
## 5 -8.879e+04
## 6 -5.801e+04
## 7 -4.737e+04
## 8 -4.421e+04
## 9 -4.351e+04
## 10 -4.336e+04
## 11 -4.335e+04
## 12 -4.335e+04
## 13 -4.335e+04
## Adaptive bandwidth (number of nearest neighbours): 24 AICc value: 86772.05
## Iteration Log-Likelihood(With bandwidth: 21 )
## =========================
## 0 -1.153e+05
## 1 -2.433e+06
## 2 -9.528e+05
## 3 -3.691e+05
## 4 -1.558e+05
## 5 -8.03e+04
## 6 -5.435e+04
## 7 -4.538e+04
## 8 -4.27e+04
## 9 -4.209e+04
## 10 -4.196e+04
## 11 -4.195e+04
## 12 -4.195e+04
## 13 -4.195e+04
## Adaptive bandwidth (number of nearest neighbours): 21 AICc value: 83978.72
## Iteration Log-Likelihood(With bandwidth: 20 )
## =========================
## 0 -1.136e+05
## 1 -2.228e+06
## 2 -8.768e+05
## 3 -3.406e+05
## 4 -1.445e+05
## 5 -7.575e+04
## 6 -5.238e+04
## 7 -4.443e+04
## 8 -4.212e+04
## 9 -4.16e+04
## 10 -4.149e+04
## 11 -4.149e+04
## 12 -4.149e+04
## 13 -4.149e+04
## Adaptive bandwidth (number of nearest neighbours): 20 AICc value: 83049.1
## Iteration Log-Likelihood(With bandwidth: 18 )
## =========================
## 0 -1.084e+05
## 1 -1.998e+06
## 2 -7.998e+05
## 3 -3.122e+05
## 4 -1.321e+05
## 5 -6.932e+04
## 6 -4.869e+04
## 7 -4.22e+04
## 8 -4.042e+04
## 9 -4.001e+04
## 10 -3.994e+04
## 11 -3.993e+04
## 12 -3.993e+04
## 13 -3.993e+04
## Adaptive bandwidth (number of nearest neighbours): 18 AICc value: 79945.88
## Iteration Log-Likelihood(With bandwidth: 18 )
## =========================
## 0 -1.084e+05
## 1 -1.998e+06
## 2 -7.998e+05
## 3 -3.122e+05
## 4 -1.321e+05
## 5 -6.932e+04
## 6 -4.869e+04
## 7 -4.22e+04
## 8 -4.042e+04
## 9 -4.001e+04
## 10 -3.994e+04
## 11 -3.993e+04
## 12 -3.993e+04
## 13 -3.993e+04
## Adaptive bandwidth (number of nearest neighbours): 18 AICc value: 79945.88gwr1 <- ggwr.basic(NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, data = ssdataspdf, bw=bw, family ="poisson", kernel = "exponential", adaptive =TRUE, dMat=DM)## Iteration Log-Likelihood
## =========================
## 0 -1.087e+05
## 1 -2.794e+06
## 2 -1.121e+06
## 3 -4.336e+05
## 4 -1.793e+05
## 5 -8.607e+04
## 6 -5.274e+04
## 7 -4.193e+04
## 8 -3.908e+04
## 9 -3.846e+04
## 10 -3.832e+04
## 11 -3.831e+04
## 12 -3.831e+04
## 13 -3.831e+04gwr1## ***********************************************************************
## * Package GWmodel *
## ***********************************************************************
## Program starts at: 2025-01-11 18:20:40
## Call:
## ggwr.basic(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV +
## POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, data = ssdataspdf,
## bw = bw, family = "poisson", kernel = "exponential", adaptive = TRUE,
## dMat = DM)
##
## Dependent (y) variable: NOR_TXCURR
## Independent variables: POP_DEN PHF_COV AHF_COV POP_POV HEALTH_DEP EDUC_DEP WORK_DEP SEC_SHOCK
## Number of data points: 123
## Used family: poisson
## ***********************************************************************
## * Results of Generalized linear Regression *
## ***********************************************************************
##
## Call:
## NULL
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.000 -0.387 -0.017 0.615 279.945
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## Intercept 8.22779520411 0.00628295755 1309.54 <0.0000000000000002 ***
## POP_DEN 0.00006017078 0.00000219293 27.44 <0.0000000000000002 ***
## PHF_COV 0.01568763771 0.00016649451 94.22 <0.0000000000000002 ***
## AHF_COV -0.00618418686 0.00003826096 -161.63 <0.0000000000000002 ***
## POP_POV -0.00000317094 0.00000005951 -53.28 <0.0000000000000002 ***
## HEALTH_DEP 0.00000701185 0.00000018534 37.83 <0.0000000000000002 ***
## EDUC_DEP -0.00003776003 0.00000031498 -119.88 <0.0000000000000002 ***
## WORK_DEP 0.00002833128 0.00000026063 108.70 <0.0000000000000002 ***
## SEC_SHOCK -0.00006747452 0.00000042839 -157.51 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 304084 on 122 degrees of freedom
## Residual deviance: 136437 on 114 degrees of freedom
## AIC: 136455
##
## Number of Fisher Scoring iterations: 5
##
##
## AICc: 136456.5
## Pseudo R-square value: 0.5513183
## ***********************************************************************
## * Results of Geographically Weighted Regression *
## ***********************************************************************
##
## *********************Model calibration information*********************
## Kernel function: exponential
## Adaptive bandwidth: 18 (number of nearest neighbours)
## Regression points: the same locations as observations are used.
## Distance metric: A distance matrix is specified for this model calibration.
##
## ************Summary of Generalized GWR coefficient estimates:**********
## Min. 1st Qu. Median 3rd Qu. Max.
## Intercept 7.5694535226 7.8182781956 8.1580282754 8.5895840674 9.2153
## POP_DEN -0.0001761390 -0.0000718419 0.0000038045 0.0001266410 0.0003
## PHF_COV -0.0293228488 0.0061415428 0.0137416914 0.0213268175 0.0291
## AHF_COV -0.0116756893 -0.0083821354 -0.0059160307 -0.0046285705 -0.0032
## POP_POV -0.0000070754 -0.0000048468 -0.0000039438 -0.0000007892 0.0000
## HEALTH_DEP -0.0000434281 -0.0000064849 0.0000054519 0.0000080731 0.0000
## EDUC_DEP -0.0000423468 -0.0000273992 -0.0000230528 -0.0000082143 0.0001
## WORK_DEP 0.0000021100 0.0000118604 0.0000317444 0.0000430481 0.0001
## SEC_SHOCK -0.0000754347 -0.0000660665 -0.0000502261 -0.0000392972 0.0000
## ************************Diagnostic information*************************
## Number of data points: 123
## GW Deviance: 75629.48
## AIC : 75692.78
## AICc : 75715.66
## Pseudo R-square value: 0.7512876
##
## ***********************************************************************
## Program stops at: 2025-01-11 18:20:40gwr$bandwidth## NULLgwr$SDF## class : SpatialPolygonsDataFrame
## features : 123
## extent : 4.978159, 9.467367, 4.270418, 7.573423 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 30
## names : Intercept, POP_DEN, PHF_COV, AHF_COV, POP_POV, HEALTH_DEP, EDUC_DEP, WORK_DEP, SEC_SHOCK, y, yhat, residual, Intercept_SE, POP_DEN_SE, PHF_COV_SE, ...
## min values : 8.2111747226557, 0.0000463682087169116, 0.0149900553967944, -0.00623666506062062, -0.00000332739638422906, 0.0000068130586121287, -0.000037262860410575, 0.0000271300062146162, -0.0000679106363728112, 0, 0.498680455500644, -4845.78772155221, 0.0443679751755873, 0.0000108111713605242, 0.000313542139110611, ...
## max values : 8.2704913487846, 0.0000615754326115069, 0.0158318110014362, -0.00604208315301457, -0.00000289996024892797, 0.00000703381128054837, -0.0000358893557158914, 0.0000288739846652074, -0.0000665217652478589, 13299, 7781.74896182731, 6429.20407371782, 0.0467326873032881, 0.0000112903996107685, 0.000342760162663926, ...gc <- round(cor(as.data.frame(gwr$SDF[,2:11]), use ="complete.obs"),2)
corrplot.mixed(gc, tl.cex = 0.8, tl.col="black")
pairs(as(gwr$SDF, "data.frame")[,2:11], pch=".")
gwr$GW.diagnostic## $gw.deviance
## 1
## 133534.9
##
## $AICc
## 1
## 133555.7
##
## $AIC
## 1
## 133554
##
## $pseudo.R2
## 1
## 0.5608618
##
## $edf
## [1] 14866.85Create spatial data frame
ssdataspdf@data$y<-gwr$SDF$y
ssdataspdf@data$yhat<-gwr$SDF$yhat
ssdataspdf@data$residual<-gwr$SDF$residual
rsd=sd(ssdataspdf@data$residual)
ssdataspdf@data$stdRes<-(ssdataspdf@data$residual)/sd(ssdataspdf@data$residual)
ssdataspdf@data$LLN=ssdataspdf@data$yhat-1.645*rsd
ssdataspdf@data$ULN=ssdataspdf@data$yhat+1.645*rsd
# Intercept
ssdataspdf@data$Intercept<-gwr$SDF$Intercept
ssdataspdf@data$est_POP_DEN<-gwr$SDF$POP_DEN
ssdataspdf@data$est_PHF_COV<-gwr$SDF$PHF_COV
ssdataspdf@data$est_AHF_COV<-gwr$SDF$AHF_COV
ssdataspdf@data$est_POP_POV<-gwr$SDF$POP_POV
ssdataspdf@data$est_HEALTH_DEP<-gwr$SDF$HEALTH_DEP
ssdataspdf@data$est_EDUC_DEP<-gwr$SDF$EDUC_DEP
ssdataspdf@data$est_WORK_DEP<-gwr$SDF$WORK_DEP
ssdataspdf@data$est_SEC_SHOCK<-gwr$SDF$SEC_SHOCK
# T-values
ssdataspdf@data$t_Intercept<-gwr$SDF$Intercept_TV
ssdataspdf@data$t_POP_DEN<-gwr$SDF$POP_DEN_TV
ssdataspdf@data$t_PHF_COV<-gwr$SDF$PHF_COV_TV
ssdataspdf@data$t_AHF_COV<-gwr$SDF$AHF_COV_TV
ssdataspdf@data$t_POP_POV<-gwr$SDF$POP_POV_TV
ssdataspdf@data$t_HEALTH_DEP<-gwr$SDF$HEALTH_DEP_TV
ssdataspdf@data$t_EDUC_DEP<-gwr$SDF$EDUC_DEP_TV
ssdataspdf@data$t_WORK_DEP<-gwr$SDF$WORK_DEP_TV
ssdataspdf@data$t_SEC_SHOCK<-gwr$SDF$SEC_SHOCK_TV
# Calculate psudo-t values
ssdataspdf@data$p_POP_DEN<-2*pt(-abs(gwr$SDF$POP_DEN_TV),df=3103)
ssdataspdf@data$p_PHF_COV<-2*pt(-abs(gwr$SDF$PHF_COV_TV),df=3103)
ssdataspdf@data$p_AHF_COV<-2*pt(-abs(gwr$SDF$AHF_COV_TV),df=3103)
ssdataspdf@data$p_POP_POV<-2*pt(-abs(gwr$SDF$POP_POV_TV),df=3103)
ssdataspdf@data$p_HEALTH_DEP<-2*pt(-abs(gwr$SDF$HEALTH_DEP_TV),df=3103)
ssdataspdf@data$p_EDUC_DEP<-2*pt(-abs(gwr$SDF$EDUC_DEP_TV),df=3103)
ssdataspdf@data$p_WORK_DEP<-2*pt(-abs(gwr$SDF$WORK_DEP_TV),df=3103)
ssdataspdf@data$p_SEC_SHOCK<-2*pt(-abs(gwr$SDF$SEC_SHOCK_TV),df=3103)
ssdataspdf$sig_POP_DEN <-ifelse(ssdataspdf@data$est_POP_DEN > 0 &
ssdataspdf@data$p_POP_DEN <= 0.05 , 'Yes', "No")
ssdataspdf$sig_PHF_COV <-ifelse(ssdataspdf@data$est_PHF_COV > 0 &
ssdataspdf@data$p_PHF_COV <= 0.05 , 'Yes', "No")
ssdataspdf$sig_AHF_COV <-ifelse(ssdataspdf@data$est_AHF_COV > 0 &
ssdataspdf@data$p_AHF_COV <= 0.05 , 'Yes', "No")
ssdataspdf$sig_POP_POV <-ifelse(ssdataspdf@data$est_POP_POV > 0 &
ssdataspdf@data$p_POP_POV <= 0.05 , 'Yes', "No")
ssdataspdf$sig_HEALTH_DEP <-ifelse(ssdataspdf@data$est_HEALTH_DEP > 0 &
ssdataspdf@data$p_HEALTH_DEP <= 0.05 , 'Yes', "No")
ssdataspdf$sig_EDUC_DEP <-ifelse(ssdataspdf@data$est_EDUC_DEP > 0 &
ssdataspdf@data$p_EDUC_DEP <= 0.05 , 'Yes', "No")
ssdataspdf$sig_WORK_DEP <-ifelse(ssdataspdf@data$est_WORK_DEP > 0 &
ssdataspdf@data$p_WORK_DEP <= 0.05 , 'Yes', "No")
ssdataspdf$sig_SEC_SHOCK <-ifelse(ssdataspdf@data$est_SEC_SHOCK > 0 &
ssdataspdf@data$p_SEC_SHOCK <= 0.05 , 'Yes', "No")#Local estimates
popden <- tm_shape(ssdataspdf) + tm_polygons(col= "est_POP_DEN", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Population Density", main.title.size = 0.75)
phfcov <- tm_shape(ssdataspdf) + tm_polygons(col= "est_PHF_COV", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "PHF Coverage", main.title.size = 0.75)
ahfcov <- tm_shape(ssdataspdf) + tm_polygons(col= "est_AHF_COV", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "ART HF Coverage", main.title.size = 0.75)
ppov <- tm_shape(ssdataspdf) + tm_polygons(col= "est_POP_POV", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Population Poverty", main.title.size = 0.75)
hdev <- tm_shape(ssdataspdf) + tm_polygons(col= "est_HEALTH_DEP", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Health Deprivation", main.title.size = 0.75)
edev <- tm_shape(ssdataspdf) + tm_polygons(col= "est_EDUC_DEP", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Education Deprivation", main.title.size = 0.75)
wdev <- tm_shape(ssdataspdf) + tm_polygons(col= "est_WORK_DEP", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Work Deprivation", main.title.size = 0.75)
sshk <- tm_shape(ssdataspdf) + tm_polygons(col= "est_SEC_SHOCK", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Security Shock", main.title.size = 0.75)
tmap_arrange(popden, phfcov, ahfcov, ppov, hdev, edev, wdev, sshk, ncol = 4)## Legend labels were too wide. Therefore, legend.text.size has been set to 0.64. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## Legend labels were too wide. Therefore, legend.text.size has been set to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## Legend labels were too wide. Therefore, legend.text.size has been set to 0.5. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## Legend labels were too wide. Therefore, legend.text.size has been set to 0.32. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## Legend labels were too wide. Therefore, legend.text.size has been set to 0.38. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## Legend labels were too wide. Therefore, legend.text.size has been set to 0.22. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## Legend labels were too wide. Therefore, legend.text.size has been set to 0.56. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## Legend labels were too wide. Therefore, legend.text.size has been set to 0.22. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
tpopden <- tm_shape(ssdataspdf) + tm_polygons(col= "t_POP_DEN", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Population Density", main.title.size = 0.75)
tphfcov <- tm_shape(ssdataspdf) + tm_polygons(col= "t_PHF_COV", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "PHF Coverage", main.title.size = 0.75)
tahfcov <- tm_shape(ssdataspdf) + tm_polygons(col= "t_AHF_COV", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "ART HF Coverage", main.title.size = 0.75)
tppov <- tm_shape(ssdataspdf) + tm_polygons(col= "t_POP_POV", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Population Poverty", main.title.size = 0.75)
thdev <- tm_shape(ssdataspdf) + tm_polygons(col= "t_HEALTH_DEP", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Health Deprivation", main.title.size = 0.75)
tedev <- tm_shape(ssdataspdf) + tm_polygons(col= "t_EDUC_DEP", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Education Deprivation", main.title.size = 0.75)
twdev <- tm_shape(ssdataspdf) + tm_polygons(col= "t_WORK_DEP", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Work Deprivation", main.title.size = 0.75)
tsshk <- tm_shape(ssdataspdf) + tm_polygons(col= "t_SEC_SHOCK", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Security Shock", main.title.size = 0.75)
tmap_arrange(tpopden, tphfcov, tahfcov, tppov, thdev, tedev, twdev, tsshk, ncol = 4)## Legend labels were too wide. Therefore, legend.text.size has been set to 0.64. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.## Legend labels were too wide. Therefore, legend.text.size has been set to 0.55. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
spopden <- tm_shape(ssdataspdf) + tm_polygons(col= "sig_POP_DEN", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Population Density", main.title.size = 0.75)
sphfcov <- tm_shape(ssdataspdf) + tm_polygons(col= "sig_PHF_COV", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "PHF Coverage", main.title.size = 0.75)
sahfcov <- tm_shape(ssdataspdf) + tm_polygons(col= "sig_AHF_COV", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "ART HF Coverage", main.title.size = 0.75)
sppov <- tm_shape(ssdataspdf) + tm_polygons(col= "sig_POP_POV", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Population Poverty", main.title.size = 0.75)
shdev <- tm_shape(ssdataspdf) + tm_polygons(col= "sig_HEALTH_DEP", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Health Deprivation", main.title.size = 0.75)
sedev <- tm_shape(ssdataspdf) + tm_polygons(col= "sig_EDUC_DEP", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Education Deprivation", main.title.size = 0.75)
swdev <- tm_shape(ssdataspdf) + tm_polygons(col= "sig_WORK_DEP", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 0.55, legend.outside = TRUE, main.title = "Work Deprivation", main.title.size = 0.75)
ssshk <- tm_shape(ssdataspdf) + tm_polygons(col= "sig_SEC_SHOCK", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Security Shock", main.title.size = 0.75)
tmap_arrange(spopden, sphfcov, sahfcov, sppov, shdev, sedev, swdev, ssshk, ncol = 4)
tm_shape(ssdataspdf) + tm_polygons(col= "stdRes", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "GWRP Std. Residuals", style = "cont")
library(writexl)## Warning: package 'writexl' was built under R version 4.2.3write_xlsx(ssdataspdf@data, "C:/Users/xyaze/Desktop/ssdataspdf.xls")library(spatialreg)## Warning: package 'spatialreg' was built under R version 4.2.3## Loading required package: Matrix##
## Attaching package: 'spatialreg'## The following objects are masked from 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOptionev <- eigenw(sscastlistw)
W <- as(sscastlistw, "CsparseMatrix")
trMatc <- trW(W, type="mult")
txlageig <- lagsarlm(NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK,
data = rschdata, listw=sscastlistw,
method="eigen", quiet=FALSE, control=list(pre_eig=ev, OrdVsign=1))##
## Spatial lag model
## Jacobian calculated using neighbourhood matrix eigenvalues
##
## rho: -0.5735677 function value: -1146.661
## rho: 0.02748166 function value: -1109.931
## rho: 0.3989506 function value: -1095.66
## rho: 0.628531 function value: -1094.018
## rho: 0.5824889 function value: -1093.742
## rho: 0.5636673 function value: -1093.725
## rho: 0.5674025 function value: -1093.724
## rho: 0.5671888 function value: -1093.724
## rho: 0.5672039 function value: -1093.724
## rho: 0.5672041 function value: -1093.724
## rho: 0.5672041 function value: -1093.724
## rho: 0.567204 function value: -1093.724
## rho: 0.567204 function value: -1093.724
## rho: 0.567204 function value: -1093.724
## rho: 0.567204 function value: -1093.724## Warning in lagsarlm(NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + : inversion of asymptotic covariance matrix failed for tol.solve = 0.000000000000000222044604925031
## reciprocal condition number = 2.55991e-18 - using numerical Hessian.## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.5672075 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.5672006 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.5672075 function value: -1093.724
## Hessian: rho: 0.5672075 function value: -1093.724
## Hessian: rho: 0.5672075 function value: -1093.724
## Hessian: rho: 0.5672075 function value: -1093.724
## Hessian: rho: 0.5672075 function value: -1093.724
## Hessian: rho: 0.5672075 function value: -1093.724
## Hessian: rho: 0.5672075 function value: -1093.724
## Hessian: rho: 0.5672075 function value: -1093.724
## Hessian: rho: 0.5672075 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724
## Hessian: rho: 0.567204 function value: -1093.724txx <- summary(txlageig, correlation=TRUE)
coef(txx)## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1341.926107743 561.250367064 2.390958 0.016804471
## POP_DEN 0.173660156 0.211841653 0.819764 0.412350670
## PHF_COV 42.583424408 16.504786307 2.580065 0.009878167
## AHF_COV -1.762526576 0.828863099 -2.126439 0.033466748
## POP_POV -0.007285443 0.004060498 -1.794224 0.072777412
## HEALTH_DEP 0.021234775 0.013142936 1.615680 0.106163586
## EDUC_DEP -0.064440087 0.021990504 -2.930360 0.003385698
## WORK_DEP 0.014471621 0.014182363 1.020396 0.307540858
## SEC_SHOCK -0.050739839 0.026472907 -1.916670 0.055279806txx##
## Call:lagsarlm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV +
## POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, data = rschdata,
## listw = sscastlistw, method = "eigen", quiet = FALSE, control = list(pre_eig = ev,
## OrdVsign = 1))
##
## Residuals:
## Min 1Q Median 3Q Max
## -3050.92 -890.83 -172.36 586.74 7774.97
##
## Type: lag
## Coefficients: (numerical Hessian approximate standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1341.9261077 561.2503671 2.3910 0.016804
## POP_DEN 0.1736602 0.2118417 0.8198 0.412351
## PHF_COV 42.5834244 16.5047863 2.5801 0.009878
## AHF_COV -1.7625266 0.8288631 -2.1264 0.033467
## POP_POV -0.0072854 0.0040605 -1.7942 0.072777
## HEALTH_DEP 0.0212348 0.0131429 1.6157 0.106164
## EDUC_DEP -0.0644401 0.0219905 -2.9304 0.003386
## WORK_DEP 0.0144716 0.0141824 1.0204 0.307541
## SEC_SHOCK -0.0507398 0.0264729 -1.9167 0.055280
##
## Rho: 0.5672, LR test value: 35.26, p-value: 0.0000000028847
## Approximate (numerical Hessian) standard error: 0.081522
## z-value: 6.9577, p-value: 0.0000000000034583
## Wald statistic: 48.41, p-value: 0.0000000000034583
##
## Log likelihood: -1093.724 for lag model
## ML residual variance (sigma squared): 2866000, (sigma: 1692.9)
## Number of observations: 123
## Number of parameters estimated: 11
## AIC: 2209.4, (AIC for lm: 2242.7)
##
## Approximate correlation of coefficients
## rho (Intercept) POP_DEN PHF_COV AHF_COV POP_POV HEALTH_DEP
## (Intercept) -0.36
## POP_DEN -0.14 -0.14
## PHF_COV 0.08 -0.70 0.12
## AHF_COV 0.27 -0.31 0.04 0.03
## POP_POV -0.05 0.07 0.31 -0.06 -0.16
## HEALTH_DEP 0.02 -0.25 -0.33 0.14 0.21 -0.70
## EDUC_DEP 0.07 -0.02 0.33 -0.01 -0.12 0.70 -0.70
## WORK_DEP -0.23 0.16 -0.15 0.02 -0.12 -0.17 -0.33
## SEC_SHOCK 0.25 -0.18 0.04 0.01 0.01 -0.22 -0.16
## EDUC_DEP WORK_DEP
## (Intercept)
## POP_DEN
## PHF_COV
## AHF_COV
## POP_POV
## HEALTH_DEP
## EDUC_DEP
## WORK_DEP -0.12
## SEC_SHOCK -0.15 -0.05txerrorlag <- errorsarlm(NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK,
data = rschdata, listw=sscastlistw)## Warning in errorsarlm(NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + : inversion of asymptotic covariance matrix failed for tol.solve = 0.000000000000000222044604925031
## reciprocal condition number = 8.18826e-18 - using numerical Hessian.txerrorlag##
## Call:
## errorsarlm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV +
## POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, data = rschdata,
## listw = sscastlistw)
## Type: error
##
## Coefficients:
## lambda (Intercept) POP_DEN PHF_COV AHF_COV
## 0.702144460 2288.142081499 0.086918521 48.297196678 -0.346993647
## POP_POV HEALTH_DEP EDUC_DEP WORK_DEP SEC_SHOCK
## -0.008551277 0.008546915 -0.049243440 0.017489395 0.001965211
##
## Log likelihood: -1097.267summary(txerrorlag)##
## Call:errorsarlm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV +
## POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, data = rschdata,
## listw = sscastlistw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3013.85 -932.99 -161.60 482.54 7734.31
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2288.1420815 689.5807901 3.3182 0.0009061
## POP_DEN 0.0869185 0.1972520 0.4406 0.6594685
## PHF_COV 48.2971967 15.5461510 3.1067 0.0018919
## AHF_COV -0.3469936 0.8161494 -0.4252 0.6707204
## POP_POV -0.0085513 0.0038551 -2.2182 0.0265426
## HEALTH_DEP 0.0085469 0.0124485 0.6866 0.4923476
## EDUC_DEP -0.0492434 0.0196173 -2.5102 0.0120660
## WORK_DEP 0.0174894 0.0175883 0.9944 0.3200407
## SEC_SHOCK 0.0019652 0.0291219 0.0675 0.9461979
##
## Lambda: 0.70214, LR test value: 28.176, p-value: 0.00000011079
## Approximate (numerical Hessian) standard error: 0.076702
## z-value: 9.1542, p-value: < 0.000000000000000222
## Wald statistic: 83.8, p-value: < 0.000000000000000222
##
## Log likelihood: -1097.267 for error model
## ML residual variance (sigma squared): 2877300, (sigma: 1696.3)
## Number of observations: 123
## Number of parameters estimated: 11
## AIC: NA (not available for weighted model), (AIC for lm: 2242.7)ssdata$pos <- predict(poisson_output)
ssdata$nbr <- predict(gnb)
ssdata$txeig <- predict(txlageig)## This method assumes the response is known - see manual pagessdata$txerr <- predict(txerrorlag)## This method assumes the response is known - see manual pagetm_shape(ssdata) + tm_polygons(col = "nbr", palette = "GnBu", title = "", border.col = "grey80", border.alpha = 0.5)
tm_shape(ssdata) + tm_polygons(col = "pos", palette = "GnBu", title = "", border.col = "grey80", border.alpha = 0.5)
tm_shape(ssdata) + tm_polygons(col = "txeig", palette = "GnBu", title = "txeig", border.col = "grey80", border.alpha = 0.5)
tm_shape(ssdata) + tm_polygons(col = "txerr", palette = "GnBu", title = "txerr", border.col = "grey80", border.alpha = 0.5)
ssdata$posres <- residuals(poisson_output)
ssdata$nbres <- residuals(gnb)
ssdata$txeigres <- residuals(txlageig)
ssdata$txerres <- residuals(txerrorlag)
tm_shape(ssdata) + tm_polygons(col = "nbres", palette = "GnBu", title = "", border.col = "grey80", border.alpha = 0.5)
tm_shape(ssdata) + tm_polygons(col = "posres", palette = "GnBu", title = "", border.col = "grey80", border.alpha = 0.5)
tm_shape(ssdata) + tm_polygons(col = "txerres", palette = "GnBu", title = "", border.col = "grey80", border.alpha = 0.5)
tm_shape(ssdata) + tm_polygons(col = "txeigres", palette = "GnBu", title = "", border.col = "grey80", border.alpha = 0.5)
Trends
trend <- read_excel("C:/Users/xyaze/Desktop/trend.xlsx", sheet = "Sheet1")
trenddata <- merge(ssmap, trend)
trenddata <- st_as_sf(trenddata)
trendy <- read_excel("C:/Users/xyaze/Desktop/trend.xlsx", sheet = "Sheet4")
trendydata <- merge(ssmap, trendy)
trendydata <- st_as_sf(trendydata)
trend1 <- read_excel("C:/Users/xyaze/Desktop/trend.xlsx", sheet = "Sheet3")
trends <- merge(ssmap, trend1)ggplot(trend1, aes(TX_CURR, fill = factor(State))) + geom_histogram() + scale_fill_viridis_d(option = "viridis", direction = -1) + theme_bw() + facet_wrap(~YEAR) + labs(x = "TX_CURR", y = "Frequency") + theme(legend.position = "bottom", legend.title = element_blank(), legend.text = element_text(size = 0.5))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 103 rows containing non-finite values (`stat_bin()`).
trend1 <- read_excel("C:/Users/xyaze/Desktop/trend.xlsx", sheet = "Sheet3")
trends <- merge(ssmap, trend1)
tm_shape(trends) + tm_polygons('TX_CURR', title = "TX_CURR", palette = "-viridis", border.col = "white", border.alpha = 0.5) + tm_facets("YEAR", free.coords = FALSE) + tm_layout(panel.labels = c('2018', '2019', '2020', '2021', '2022'), legend.show = FALSE)
` ##